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1.0 PROGRAM PLAN OUTLINE AND NARRATIVE 

1.1 Introduction 

This Annual Report focuses on the effort that has been completed 
during the second year of the technical effort. The total project Is now 
expected to last a total of five years. All elements of the technical 
tasks to be accomplished have now been defined. The new effort in the 
third year Is the Initiation of programming of the advanced methods 
formulation; the approximate methods effort will not begin until the fourth 
year of technical effort. 

The project Is very much a team effort with significant contributions 
coming from several task managers: Dr. O.H. Burnside - Verlflcatlon/Valldatlon 
Coordination; Dr. Y.-T. Hu - NESSUS/FPI Development; Drs. J. Nagtegaal, S. 
Nakazawa and Mr. J. Diaz - PFEM Development; Or. K.R. Rajagopal - Verification 
Studies; Dr. P. Fink - NESSUS/EXPERT Development; Prof. P. Wlrschlng - 
Advanced Simulation Methods; and Prof. S. Atlurl - Hybrid FEM Development and 
Level III Modeling. The SwRI Program Manager acknowledges the critical 
contributions from each of these Individuals. 

The remainder of this Section outlines the elements of the technical 
approach being taken In PSAM. Section 2.0 summarizes the technical 
accomplishments of the second year of the project, supported by various 
appendices. Section 3.0 presents a brief outline of some of the current 
efforts. 

1.2 Probabilistic Finite Element Methods (PFEM) Plan 

The developed methods of analysis are to treat linear problems as well 
as those with nonlinear material and geometric response. Stochastic 
modeling of loads (e.g., centrifugal, thermal, pressure), geometry, and 
material behavior are being modeled with three levels of approximation, 
relative to accuracy and confidence. Level I analyses treat randomness as 
being spatially homogeneous (e.g., each part has a different modulus, yield 
stress, thermal load, etc.). Level II analyses treat random variables as 
random fields (e.g., modulus variability Is different in the bore of a 
disk, versus the rim of the disk; pressure uncertainty Is different at the 
root of an airfoil versus the tip of the blade). Level III stochastic 
modeling is to be able to reflect uncertainty between variables In the 
governing equations (e.g., strain agrees with displacement gradients only 
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In stochastic, not deterministic, terms; stress Is related to strain 
through stochastic relations). 

Two methods of probabilistic modeling are Included In the various 
analysis methods. The first of these Is the Fast Probability Integration 
(FPI) method. The FPI method Is adopted from the field of structural 
reliability as a way of predicting the probability that a response variable 
(e.g., stress, frequency) will exceed some allowable. The method Is based on 
establishing the approximate sensitivity of a response variable to the 
stochastic variables, and then processing these sensitivities by the FPI 
algorithm to establish the cumulative response distributions for the 
variables. The second method Is direct simulation using an enhanced Mon te 
Carlo method. Both probabilistic prediction methods will make use of the same 
structural s en sitivity data ba se, whic h Is generated by NESSUS. Confidence 
levels will be estimated for the response variable distributions that are 
calculated. A composite load spectrum analysis procedure will be Included. 

The PFEM Is a direct adaptation of standard finite element methodology to 
the needs of PSAM. The finite element code (NESSUS/FEM) Is to Include 
plate and shell elements based Initially on the displacement method of 
formulation, and on linear equations of motion and material behavior. 

Hybrid plate and shell elements are to be Included, as well as nonlinear 
geometric and material behavior. The NESSUS/FEM code will Include a 
variety of standard finite elements for structural modeling. The 
NESSUS/FEM program will allow for nonlinear el astoplastl c/creep modeling, 
and for geometric nonlinear problems of finite displacement, rotation, and 
strain. 

In addition, an enhanced shell/plate element formulation will be 
developed. This enhanced formulation will be a quasi -continuum element 
that provides for surface data Input and nodal stress recovery, consistent 
with the requirements of the NASA SOW. The enhanced element Is a 
displacement formulation, developed from the Hu-Washlzu variational 
formulation. Stresses, strains, and displacements will be Interpolated 
Independently. In order to reduce the formulation to a displacement- like 
formulation, the stress and strain fields are discontinuous between the 
elements. Displacements will be Interpolated on a nodal basis, with nodes 
selected at the surfaces of the shell/plate element. The element will 
satisfy all constant stress modes and will provide full rigid body mode 
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capability (l.e., has correct rank). The eight noded element will provide 
for surface pressure load definition, as well as for nodal stress, strain 
recovery. 

The hybrid element formulation Is also based on the Hu-Washlzu 
variational statement. Thus, It will have stress modes that are defined 
Independent of the displacement modes. The element will be a sixteen node 
shell/plate element with surface loading and nodal stress, strain recovery 
capability. Special interpolation capability for severe thermal gradients Is 
planned In both the enhanced and hybrid shell/plate formulations. 

Material response Is to Include the range from elastic to 
thermoviscoplastic. The material model will be based on theoretical 
development for a random relationship between stress and strain for the 
general class of thermomechanlcal response problems. The material modeling 
considerations will allow for a full. Level III Interpretation of 
stress/straln stochastlclty. The model will be based on the assumption 
that each material has its own stochastic response over the full range of 
loading history. Thus, we rule out as a mathematical construct, the notion 
of Incremental stochastlcity. The theoretical material modeling development 
will admit Implementation of endochronlc or "thermoviscoplastic considerations. 

The NESSUS code Is modular for adaptation to the General Purpose 
Structural Analysis (GPSA) framework. The modules Include NESSUS/FEM, 
NESSUS/PAAM, NESSUS/BEM, NESSUS/FPI, NESSUS/PRE, NESSUS/EXPERT, and 
others as needed. Interfaces between these modules will be clearly 
defined. 

The approach to validation Is to perform validation and verification 
studies on the new element and formulation capabilities as they become 
available. This will also provide for direct comparisons between the 
various solution capabilities. The NESSUS code Is to be continually 
validated through Its application to well-defined problems with known 
probabilistic responses In order to demonstrate the full and reliable 
capability of the code. 

It has been found to be very important that the verification study 
Include a wide range of simple structural models that exercise the various 
options of the NESSUS code. These verification problems, being run by SwRI 
and Rocketdyne, serve to provide further confidence In the code, to develop 
rule bases for NESSUS/EXPERT, and demonstrate the utility of the code. 


4 


The NESSUS code will be verified by Its application to four selected 
space propulsion system hardware Items. These will Include the turbine 
blade, transfer duct, LOX post, and the high pressure oxidizer duct. 
Experimental data to support the analyses will be compiled and 
statistically modeled. The four verification problems are outlined In 
Appendix A. ^ ...... iS .. .... - ... 

1.3 Probabilistic Approximate Analysis Methods (PAAM) Plan 

The PAAM code will be established by SwRI In consultation with 
Rocketdyne staff. The purpose of the PAAM code Is to provide a mechanics 
of materials approach to the probabilistic modeling of plate and shell type 
structures. The approach to be taken by SwRI will be to: 

1. Identify simplified plate/shell problems representative of plate 
and shell regions within the four selected propulsion system 
components. 

2. Identify plate and shell type analytical solutions that 
correspond best to the physical problems Identified In 1., above. 

3. Modify the analytical solutions to account. In a suitable and 
approximate manner, the loading, material response, and structural 
response features required for the four component problems. 

4. Program NESSUS/PAAM to Include a library of these solutions and 
approximation methods for loading, material response, and structural 
response. 

1.4 Probabilistic Advanced Methods (PAdvAM) Plan 

The basis of the Probabilistic Advanced Analysis Methods Is the boundary 
element method, specifically the BEST3D code previously developed under NASA 
HOST funding. SwRI has further developed this code and proposes, to modify It 
In a manner suitable for Inclusion In the PSAM analysis library as NESSUS/BEM. 

The boundary element method (BEM) contrasts, for the linear problem, 
with the finite element method (FEM) by the fact that the governing 
equations are written at the boundary of the body only. The so-called 
boundary integral equation (BIE) governs the relationship between tractions 
and displacements at the surface of the body. The only geometric 
description of the body that Is required is the surface of the body. For 
the thermoelastic problem with variable material properties and problems of 
linear vibration, It Is also possible to reduce the continuum problem to a 
boundary formulation. For problems with geometric or material nonlinearities, 
and for transient dynamic problems, a volume modeling Is generally required. 


The perturbation algorithm will be developed for NESSUS/8EM. For 
those problems with no volume Integrations, all perturbations will be In 
terms of surface data. Geometry, for example, will be perturbed through 
explicit differentiation of the boundary modeling shape functions. The 
perturbations will maintain continuity of boundary shape by moving the 
boundary node locations. Perturbations of the mass matrix for vibration 
analysis will be similarly modeled. Level I material perturbations will be 
explicitly accounted for. 

Level II and Level III material perturbations will be examined using 
one of two possible approaches. The most direct Is to handle these through 
volume Integrals (discussed below). The most Interesting Is to develop 
boundary models that can Interpolate volumetric changes. In terms of 
perturbed boundary data. The latter approach Is favored and will be the 
first to be explored. Explicit differentiation or differencing of boundary 
data will be used In order to avoid an Iterative solution algorithm. 

Volume Integration methods will be especially developed for NESSUS/BEM 
to take advantage of plate/shell type of behavior. Simplifying 
Interpolation assumptions will be made to reduce the need for significant 
numbers of volume discretizations In the through-thickness direction. It 
will be assumed that deviations In strain behavior In this direction from 
the linear solution are not excessive. 

The first year (FY87) will establish the linear static and dynamic 
thermoelastic modeling capability for NESSUS/BEM. The second year (FY88) 
will focus on the establishment of the essential nonlinear modeling 
capability, but without the full Level III thermoviscoplastic modeling and 
random transient loading. The third and final year (FY89) will release the 
full nonlinear capability. 

The stochastic basis of a variational model of structural response 
will be established by GIT researchers under the direction of Professor 
Satya Atluri. The stochastic variational statement will be used to 
demonstrate the formulation basis of the Level I, II NESSUS PFEM models. 
Further, It will be used to establish the Level III formulation for 
adoption Into NESSUS. It Is expected that the Level III model will be 
based on the use of correlation model matrices linking the strains and the 
displacement gradients, and another linking stress to the material 
constitutive behavior. 
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2.0 TECHNICAL PROGRESS SUMMARY 

2.1 Task I: PFEM 

2.1.1 NESSUS/FEM Development 

2. 1.1.1 Status at End FY85 

The finite element analysis module NESSUS/FEM has evolved 
from the MHOST code, developed by MARC for Pratt and Whitney Aircraft Company 
under NASA contract NAS3-23697. A review of the capabilities of MHOST by SwRI 
Indicated the need for enhancements to provide additional features relevant 
to the analysis of reusable space propulsion system components. This 
enhanced version of the MHOST code was delivered to SwRI In August 1985 as 
NESSUS 0.1, and was the latest version of the code shipped from MARC prior 
to the end of FY85. 

The major enhancements provided with NESSUS 0,1 Included: 

A. Element library and problem modeling features 

o Addition of a two-noded Timoshenko beam element 
o Rotational Inertia terms In consistent mass matrices 
o Grounded springs of prescribed stiffness 
o More convenient definition of time-histories for pulse 
loading 

B. Algorithmic enhancements 

o Displacement method option for linear, elastostatics 
o Power shift option for eigenvalue extraction 

C. Analysis capabilities for linear systems 

o Transient dynamics using mode superposition 
o Harmonic loading and base excitation 
o Random vibration (PSD) analysis 

By the end of FY85, the basic formulation for probabilistic finite 
element analysis as implemented in NESSUS had been developed and 
demonstrated on a few sample problems. The original approach relied on a 
Taylor series expansion of the stochastic problem about a deterministic 
solution. This approach did not appear to be practical for the large 
systems of finite element equations parameterized by many random variables 
that are needed for realistic SSME applications. An alternative approach 
was developed, based on an iterative perturbation analysis method that uses 
the factorized stiffness of the unperturbed system as the iteration 
preconditioner for obtaining the solution to the perturbed problem. This 
approach eliminates the need to compute, store and manipulate explicit partial 
derivatives of the element matrices and force vector, which not only reduces 


memory usage considerably, but also greatly simplifies the coding and 
validation tasks. A similar approach for the solution of the perturbed 
symmetric elgenproblem was developed by Professor Juan Slmo, at the Applied 
Mechanics Division, Stanford University, for Implementation In NESSUS/FEM. 

The efficient treatment of correlated random variable fields was 
Identified early on In the PSAM effort as a major practical Issue, since 
many SSME applications Involve random variables that are correlated to some 
degree. Examples of this Include random pressure and temperature fields 
defined on the surface of a turbine blade, or the thickness of the walls 
and liners in the transfer ducts or nozzle of a rocket engine. The 
strategy adopted In the NESSUS code relies on a variable transformation 
Into the elgencoordl nates of the covariance matrix defining the random 
field. The transformed variables can be shown to be uncorrelated and may 
therefore be manipulated as such In the NESSUS/FEM and NESSUS/FPI modules. 
The computation of the transformed variables may be carried out prior to 
the finite element analysis of the model and may, therefore, be regarded as 
pre-processing operation. 

Several aspects of the proposed formulation were demonstrated on an 
ad-hoc basis before the end of FY85. The feasibility of the Iterative 
perturbation algorithm for elastostatics was demonstrated in April 1985 
with a problem Involving a clamped square plate under uniform pressure 
loading, using a 10 x 10 mesh of shell elements and subjected to thickness 
variations along one edge. The numerical manipulations proposed for 
handling correlated data were demonstrated also In April 1985 with a 
problem Involving a scalar random variable field defined on a 10 x 10 grid 
with varying strength of correlation. Finally, all Ingredients for the 
proposed formulation were combined in a demonstration problem using a 
simplified model of a curved turbine blade discretized with 48 shell 
elements, and having random pressure and temperature fields with partial 
correlation, random uniform thickness, and random stiffness at the root. 
This exercise was completed In May 1985. Although the formulation for the 
iterative solution of the perturbed symmetric elgenproblem was essentially 
complete by the end of FY85, no demonstration problems using this approach 
were available at the time. 
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2. 1.1.2 Database Development 

The perturbation database (Fig. 2.1) provides an external 
record of the perturbation data obtained during execution of the NESSUS/FEM 
module. In a typical NESSUS/FEM execution, a number of perturbed solutions 
about a deterministic state are computed with the use of appropriate numerical 
algorithms. The results for both the unperturbed and all perturbed systems 
are added to the perturbation database as soon as each converged solution 
becomes available. The Information stored In the database may then be 
accessed by the NESSUS/FPI module, in order to extract the data required for 
the computation of a system reliability estimate or to obtain distribution 
curves for relevant response variables. The perturbation database Is problem- 
specific, and was designed to centralize all the Information pertinent to the 
analysis of a given model, even If It Is obtained In the course of multiple 
NESSUS/FEM executions. Future releases of NESSUS/FEM will allow full use of 
these capabilities. 

The perturbation database resides In a binary 
(unformatted) direct-access file, and may be accessed using standard FORTRAN 
I/O facilities. The database is structured as a two-way ordered linked list 
(Fig. 2.2), allowing for quick and efficient traversal In search of a specific 
entry. This type of data structure allows for the Insertion and deletion of 
Individual entries anywhere in the list without violating the original 
ordering convention (Fig. 2.3). It Is therefore possible to enrich the 
existing database with Information obtained In multiple executions of 
NESSUS/FEM without having to regenerate data obtained in previous runs. 

The organization of a typical database constructed by 
NESSUS/FEM is outlined in Fig. 2.4. The entry point Is a single PROBLEM 
HEADER RECORD, always occupying the first physical record In the file. This 
record contains sizing Information pertinent to the problem, together with 
pointers to two distinct ordered linked lists. One list contains the load 
Incrementation history, with the Individual perturbation data sets nested 
Inside each increment. The second list contains the eigenvalue and 
eigenvector data, this time with the individual elgenpalrs nested Inside each 
perturbation data set. Both lists consist In a series of INCREMENTAL or 
EIGENPAIR DATA HEADER RECORDS, forming two two-way ordered linked lists, shown 
In Fig. 2.4 as extending downward and upward from the single entry point. 

These headers in turn contain pointers to the actual DATA RECORDS, containing 



Fig. 2.1 The NESSUS Perturbation Database 
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Eigenpair 2 
Perturbation 1 


Eigenpair i 
Perturbation l 


Eigenpair 2 
Perturbation 0 


Eigenpair 1 
Perturbation 0 


Problem Header 
(Entry Point) 


Increment 0 
Perturbation 0 


Increment 0 
Perturbation 1 


Increment 0 
Perturbation 2 


Increment 1 
Perturbation 0 


Increment 1 
Perturbation 1 


Fig. 2.4 Data Structure of the NESSUS Perturbation Database 
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Information on the type of perturbation and the perturbed system response. A 
null pointer Is used to flag the unavailability of data, which may result from 
the lack of a converged solution. Null pointers are also used to terminate 
both the Incremental and elgenpalr data lists. 

The present Implementation of the perturbation database 
provides easy and efficient data retrieval using standard algorithms for 
manipulating ordered linked lists. Insertion and deletion of Individual 
entries can be accomplished locally without the need for moving large blocks 
of data. The Internal data structure was designed with the flexibility to 
accommodate additional capabilities planned for future releases of NESSUS/FEM 
with minimal adjustments to the software already In place. The use of binary 
(unformatted) files provides compact storage for the potentially massive 
amounts of data required for the analysis of realistic problems. For small 
problems, a simple FORTRAN utility Is available to provide translation of the 
database Into formatted (printable) form. This can be quite useful for 
debugging codes written to access the database, or for moving small databases 
across different computer systems. The Internal data structure of the 
perturbation database is well documented in a report which can be used as a 
guide for the development of new codes requlrihg access to existing 
databases. Finally, It should be noted that the information contained In the 
perturbation database may be useful for applications other than probabilistic 
structural analysis, such as the Investigation of the sensitivity of the 
response to several design parameters. 

2.1. 1.3 NESSUS/PRE Module Development 

The NESSUS/PRE module (Fig. 2.5) Is a pre-processor used 
for the preparation of the statistical data needed to perform probabilistic 
finite element analysis with NESSUS/FEM. NESSUS/PRE allows the user to 
describe a spatial domain defined by a set of discrete points, typically 
corresponding to the nodal points of an existing finite element mesh. One or 
more random variable fields may then be specified over this spatial domain by 
defining the mean value and standard deviation of the field variables at each 
at each point, together with the appropriate form of correlation. Each random 
variable field may be modeled as uncorrelated, fully correlated or partially 
correlated. The current version of NESSUS/PRE limits the treatment of 
partially correlated fields to fields of Gaussian variables with equal 
correlation strength in all directions (isotropic correlation). 
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Fig. 2.5 The NESSUS/PRE Module 
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Each random field Is treated according to the form of 
correlation specified for It. Uncorrelated random fields are automatically 
decomposed Into a set of uncorrelated vectors, each corresponding to a unit 
variation at a given degree-of-freedom. Fully correlated fields are 
automatically converted to a single vector, corresponding to a scaling of the 
random field by one standard deviation. For a partially correlated random 
field, the preprocessing operation In NESSUS/PRE Is considerably more 
complex. This will Involve the construction of the variance-covariance matrix 
for the field, followed by the spectral decomposition of this matrix. The 
field data Is then transformed Into the elgencoordl nates of the covariance 
matrix, yielding a set of mutually uncorrelated random vectors which contain 
all the Information present In the original correlated field. The theoretical 
details of this procedure are given In Section 5.3.1 of the PSAM First Annual 
Report. The spectral decomposition of the covariance matrix is performed 
conservatively, using Jacobian Iteration to solve simultaneously for all 
eigenvalues and eigenvectors of the matrix. If the correlation Is strong, the 
uncertainty In the data Is dominated by just a few of the highest eigenvalues 
of the matrix. Hence, the user Is given the option to simplify the problem by 
truncating the spectrum to a prescribed tolerance, retaining only the most 
significant eigenvalues for the analysis. This strategy can produce a very 
significant reduction in the amount of computation required for the analysis, 
especially In problems involving a large number of random variables. In all 
cases, the output from NESSUS/PRE will consist of a set of uncorrelated random 
vectors written to an external formatted data file. This file will contain 
the random variable definitions for NESSUS/FEM, and may be Included In the 
Input deck to the finite element module without further modification. 

The present Implementation of NESSUS/PRE allows the 
specification of random fields Involving: 

1. nodal coordinate data 

2. nodal shell thickness 

3. nodal shell or beam normals 

4. thickness of plane stress elements 

5. modulus of elasticity 

6. Poisson's ratio 
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7. thermal expansion coefficient 

8. material density 

9. rotational speed 

10. nodal force vectors 

11. element pressures and edge tractions 

12. nodal temperatures 

13. elastic beam section properties 

14. base spring stiffnesses 

15. orientation of anisotropy axes 

Additional types of random variables will be Included In 
future releases of NESSUS/PRE as required by the enhanced capabilities of 
NESSUS/FEM. 

2. 1.1.4 Code Structure 

The NESSUS/FEM module (Fig. 2.6) provides finite element 
modeling and analysis capabilities for probabilistic structural analysis 
problems. The finite element code Is structured as a set of six major driver 
routines, reflecting the types of analysis currently available. These 
Include: 


1. A static analysis driver for the solution of linear and 
nonlinear problems In either a purely Iterative manner or In 
incremental-iterative fashion. 

2. A bifurcation buckling driver, used for stability analysis of 
linearized structural systems. 

3. A modal extraction driver for the determination of the 
undamped natural frequencies and mode shapes for vibrating 
structures. 

4. A mode superposition driver for the analysis of steady state 
or transient linear vibration problems In the time domain. 

5. A random vibration driver for the analysis of problems 
Involving stationary random excitation by Integration In the 
frequency domain. 

6. A direct time integration driver using the Newmark-b method 
for the solution of linear and nonlinear transient dynamics 
problems. 
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A new facility has been added at the topmost level of 
NESSUS/FEM to allow conditional transfer of control between driver routines. 
This feature allows the performance of more sophisticated types of analysis, 
which are useful for the realistic modeling of typical SSME components. A 
typical application might Involve the static analysis of a spinning turbine 
blade, with centrifugal loading applied over a number of load Increments. 

At a prescribed Increment number, control of the execution may be 
transferred to the modal extraction driver. In order to determine the 
vibration characteristics of the blade. Including centrifugal mass and 
stress stiffening effects due to the Initial stresses obtained In the 
static analysis. These features were first available In Version 1.2 of the 
code. 

Significant efficiency Improvements were achieved by 
replacing the old band and frontal equation solvers with a newly developed 
profile solver. The new solver, available In Version 1.3, not only provided 
Increased speed In the factorization and back substitution phases of the 
analysis, but also resulted In a substantial reduction of memory requirements 
for medium to large problems. This allowed the In-core solution of large 
turbine blade models using 8-noded bricks. Selected performance results for 
the new solver are summarized In Fig. 2.7 - Fig. 2,9. These numbers were 
obtained on the PRIME 9955 at MARC, with the memory requirements expressed In 
single precision (32 bit) words. 

The extraction of eigenvalues for both linear dynamics 
and buckling problems is performed using the subspace Iteration method. 
Multiple power shifts may be used to extract modes within prescribed frequency 
bounds. This technique Is particularly useful In the analysis of structures 
containing rigid-body modes. The eigenvalue analysis subsystem Is very 
similar to the one available In NESSUS 1.0, having been modified to 
accommodate profile storage for the stiffness and mass matrices, together 
with other minor efficiency Improvements. 

A full library of modern algorithms for nonlinear 
analysis Is available In NESSUS/FEM. Both full Newton and modified Newton 
iteration algorithms have been available since Version 0.1. Newly Implemented 
algorithms for nonlinear analysis include the line search algorithm, 

Davldon rank-one secant Newton update and Inverse BFGS rank-two update. 

These algorithms have been available In Version 1.0 and up. Variations of 
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Fig. 2.7 SSME HPFTP Blade Model with 1025 Brick Elements and 1575 Nodes 


Profile Solver Band Solver 

Memory Requirement 3483811 4459647 

Solution Time 1466.594 sec N/A* 

* Too Large to Run on PRIME 9955 at MARC 
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Fig. 2.8 Buckling Analysis of a Cylinder with 160 Shell Elements and 176 
Nodes 


Memory Requirement 

(a) Static 

(b) Eigenvalue 

Solution Time 

(a) Static 

(b) Eigenvalue 


Profile Solver Band Solver 


498873 596703 

1044773 1340375 


62.057 sec 171.430 sec 

66.788 sec 187.551 sec 
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Fig. 2.9 Modal Analysis of a Composite Laminate Fan Blade with 240 Shell 
Elements and 279 Nodes 


Memory Requirement 
Solution Time 


Profile Solver 

824341 
24.794 sec 


Band Solver 

1054259 
93.764 sec 
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many of these algorithms have since been applied to the computation of the 
solution to the perturbed elastostatlc problem, as discussed In Section 
2. 1.1.6 of this report. 

2. 1.1.5 Element Technology 

The element library currently available In NESSUS 
consists of six Isoparametric, numerically Integrated element types (see Table 
2.1). Geometric quantities and material properties are defined at the nodes, 
and Interpolated Into the Interior of each element using the appropriate shape 
functions. A nodal projection and smoothing algorithm Is used to allow the 
reporting of strains and stresses on a nodal basis. 

Continuum- type problems may be modeled using bilinear 
four-node quadrilaterals for plane stress, plane strain or axlsymmetrlc 
situations, or trl linear eight-node bricks for three-dimensional problems. 

All B-matrlx routines for continuum elements allow full, reduced, trapezoidal 
and selective Integration. Selective Integration Is Implemented using the 
B-bar approach, and has been designed to facilitate the Implementation and 
testing of different Integration weighting schemes. The performance of 
these elements has recently been Improved with the adoption of a strain 
filtering scheme based on a local element orthogonal coordinate system 
constructed by polar decomposition of the Jacobian matrix for the 
Isoparametric mapping. This technique enhanced the behavior of the element 
In situations Involving distorted elements. 

The shell element currently available In NESSUS Is a 
four-node Isoparametric formulation derived from the Relssner-Mindlln plate 
and shell theory. Bilinear Interpolations are used for the coordinates, 
displacements and rotations. The element Is selectively Integrated, and 
stabilized by hourglass control on the transverse shear terms. An In-plane 
twist term Is Included to avoid the "drilling mode" singularity on a flat 
assembly of elements. This element may be used to model thick shell 
problems, with significant transverse shear deformation, and retains 
acceptable accuracy when used to model thin shell structures. 

A two-node linear Isoparametric beam element is also 
available, based on Timoshenko beam theory. Linear Interpolations are used 
for the cross-section, displacements and rotations. Reduced one-point 
Integration Is used for economy, since this will yield a rank-sufficient 
stiffness matrix for the element. Since the cross-sectional properties for 
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Table 2.1 

Summary of the NESSUS Element Library 


P.STRS P.STRN AXSYM BRICK SHELL T. BEAM 


ITYPE 

NELCRD 

NELNFR 

NELNOD 

NELSTR 

NELCHR 

NELINT 

NELLV 

NELLAY 

HD I 

NSHEAR 

JLAW 



ITYPE Element type number. 

NELCRD Humber of coordinate data per node. 

NELNFR Number of degrees -of- freedom per node. 

NELNOD Number of nodes per element. 

NELSTR Number of stress and strain components per node. 

NELCHR Number of material property data for the element. 

NELINT Number of 'full' integration points per element. 

NELLV Number of distributed load types per element. 

NELLAY Number of layers of Integration through the thickness of the 

shell element. 

NDI Number of direct stress components. 

NSHEAR Number of shear stress components. 

JLAW Type of the constitutive equation. 
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this element are defined In pre- Integrated form. Its use Is restricted to 
linear elastic problems. 

In recent PSAM meetings a strong desire has been 
expressed for the development of advanced element technology needed to address 
specific SSME applications In an effective manner. Many of these applications 
Involve localized effects which cannot be captured using the classical plate 
theories. Examples Include strong curvature, strong thickness variations 
or localized mechanical or thermal loading. In principle, continuum theory 
will always be able to model the proper solution. However, regular 
continuum elements lack the appropriate deformation modes to model 
shell-like behavior In a satisfactory way. Recent developments In element 
formulation suggest that It may be possible to construct continuum elements 
with enhanced bending behavior that would perform well when degenerated In 
one direction to form a shell-like element. The development of such an 
element was proposed by MARC for Implementation In the NESSUS code. 

2. 1.1.6 Solution Strategy and Algorithms 

The use of FPI methods In probabilistic finite element 
analysis Involves the repeated computation of the structural response for 
small perturbations of the random parameters about a given deterministic 
state. Probabilistic models of realistic structural systems can be quite 
complex, requiring the analysis of large finite element models parameterized 
by many random variables. The computational effort expended In the generation 
of perturbed solutions for these models vastly exceeds that required for all 
other phases In the analysis. Hence, the ability to efficiently compute 
the response of the perturbed system Is crucial to the viability of the 
method . 

For linear elastostatlcs, the basic perturbation problem 
may be expressed as follows. Given the solution to the unperturbed set of 
finite element equations 
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• where 
K » K + dK 


u = u + du 
f * f + df 


( 3 ) 


Substitution of these definitions Into the equation for the 
perturbed problem will yield 

(K + dK)(u + du) » (f + df) 

K du » (f + df) - (K + dK) u - dK du 

A A 

K du * f - K u - dK du 

Several methods have been proposed for the solution of 
the problem In this form. A first-order perturbation method may be obtained 
by neglecting the last term (second-order), and solving for a first-order 
approximation to du. This approximation can be shown to correspond to the 
first term In the Taylor series expansion for du. Higher-order 

perturbation methods are obtained by carrying along additional terms In the 
Taylor series expansion. 

The perturbation strategy adopted In NESSUS/FEM Is based 
on the recovery of the higher-order terms by an Iterative process. A suitable 
algorithm Is provided by the recursion form 

K du*"* 1 ) - f - K „("> (7) 

i*"* 1 ) • u<"> * di< n+1 > 

( 8 ) 

This process Is equivalent to a modified Newton 
Iteration, and can be shown to satisfy the appropriate consistency 
condition. The stability of the algorithm will be discussed In Section 
2. 1.1. 7 In some detail. 


( 5 ) 

( 6 ) 
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Many advantages may be derived from adopting this 
strategy. It should be noted that all computations with the perturbed system 
can be performed at the element level, and only the resulting element Internal 
force vectors need be assembled Into a global vector. This residual force 
vector must be computed at every Iteration step, and provides a direct measure 
of the quality of the approximation that Is used to establish convergence. 
Furthermore, this approach eliminates the need to compute and store 
explicit partial derivatives of the element stiffness and load vector, or 
any assembled form of these quantities. This not only significantly 
reduces data storage requirements, but also greatly simplifies the coding 
and validation tasks. The perturbation of geometry and material data Is 
made Independent of the element formulation adopted, which allows simple 
extension of the method to newer element technologies. The overall 
efficiency of the method can easily be justified on the basis of well-known 
operation count statistics for large finite element problems. 

Additional efficiency Improvements are obtainable from 
recasting the perturbation problem as an Iterative process. This allows the 
Implementation of a number of convergence acceleration methods for 
Iterative problems, such as the line search algorithm and quasi-Newton 
iteration schemes. In particular, significant performance Improvements 
have been demonstrated with the use of either Davldon rank-one secant 
Newton update or Inverse BFGS rank -two update applied to the perturbed 
elastostatlc problem. The present Implementation of the line search 
algorithm does not appear to be very cost-efficient for linear elastostatlc 
problems. This Is, In part, due to the fact that It has been Implemented 
as a truly nonlinear line search, since the final goals of the PSAM project 
call for the extension of the perturbation algorithms to nonlinear 
problems. 

A similar perturbation algorithm for the symmetric 
elgenproblem with iterative Improvement also has been developed and 
Incorporated in NESSUS/FEM. This algorithm differs from earlier elgenproblem 
perturbation methods by the fact that it has been developed from the start 
with the intent to tackle realistic structural vibration problems. The 
problem of properly splitting eigenvalue clusters In the spectrum of the 
unperturbed problem was Identified early on in the algorithm development. A 
solution was developed, Involving a reduced elgenproblem with the dimension of 
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the multiplicity of eigenvalues In the cluster.’ Similar formulations 
typically Involve the solution of a larger eigenvalue problem, with dimensions 
of at least the number of perturbed eigenvalues. An error estimate to account 
for the effect of a truncated modal representation Is also provided. The 
final form of the algorithm Is Independent of the method used to obtain the 
starting elgenpalrs, and can be used with any of the modern algorithms for 
the solution of large symmetric positive-definite elgenproblems. 

These perturbation analysis algorithms have been 
available In NESSUS/FEM since Version 1.0 and have been successfully applied 
to a broad class of linear problems. It Is expected that much of the code 
developed for the perturbation of linear elastostatlcs problems will be able 
to handle weakly nonlinear situations with only minor modifications. 

2. 1.1. 7 Stability Considerations 

The Iterative perturbation analysis algorithms available 
In NESSUS have been successfully applied to a broad range of structural 
problems over the past year. The experience acquired In this testing and 
validation phase also Identified a class of problems for which the Iterative 
process was observed to become unstable with seemingly small values of the 
perturbation parameter. The problem was first encountered In the analysis 
of validation problem #2, described in Section 4.0 of Volume III of the PSAM 
First Year Progress Report. This problem Involved the analysis of a thin 
cantilever beam using bilinear Relssner-Mindlln shell elements (NESSUS 
element 75) under bending. The stiffness equations for this problem are 
poorly conditioned, as a result of the enforcement of the transverse shear 
constraint In the thin limit of the Relssner-Mindlln theory. The problem 
was observed to be particularly sensitive to geometry perturbations 
Involving changes in element length, which often resulted In loss of 
stability of the Iterative algorithm even for small elongations of the 
mesh. 

A detailed Investigation Into the nature of the problem 
was undertaken at MARC, and the major findings are summarized In Appendix B to 
this report. The Investigation concentrated on the analysis of a simpler 
model problem, Involving a mesh of linearly Interpolated Timoshenko beam 
elements. This Is the one-dimensional analog of the Relssner-Mindlln plate 
problem, and exhibits pathological behavior Identical to that observed In 
the plate problem, while offering a much simpler formulation and far more 
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amenable to a detailed analysis. Early on In the course of the 
Investigation the problem was found to be governed by the characteristics 
of the assembled stiffness equations. A von Neumann stability analysis of 
the assembled equations at a typical Internal node for the model problem 
was performed, which provided closed-form expressions for the stability 
limit In the case of uniform mesh elongation. These stability limits were 
found to accurately predict loss of stability for numerical experiments 
Involving both beam and shell element discretizations of the model problem. 
These results can be used to estimate the stability limits for more general 
beam and plate problems, providing an upper bond for the size of the 
perturbation parameter that will preserve the stability of the algorithm. 

In general, stability will present a concern for the 
analysis of any problem Involving some form of Implicit constraint equations 
in the underlying theory. Stability problems will typically arise whenever 
geometry perturbations affecting these constraint equations are Imposed on 
the unperturbed problem. Such problems Include the analysis of 

o Thin plates and shells allowing shear deformation 
o Incompressible elasticity, e.g., rubber-like materials 
o Strongly anisotropic materials 
o Devlatorlc rate- Independent plasticity 
o Incompressible Stokes and Navler-Stokes flow 

Several of these problems are relevant to SSME 
applications. It must be noted that alternative formulations based on a 
Taylor series expansion about the unperturbed system are not Immune to the 
problem. This may be concluded by noting that the speed of convergence of the 
Iterative algorithm is closely related to the error associated with the 
truncation of the Taylor series. The analysis of the general problem Is 
complicated by the fact that stability Is often governed by the lowest 
deformation modes present in the assembled stiffness equations. Thus, the 
development of general closed-form results for unstructured, multi-dimensional 
meshes subjected to non-uniform distortion does not appear to be practical. 
However, the insight obtained from the analysis of simple model problems 
can be used to develop "smart" algorithms capable of adaptively adjusting 
the perturbation size in order to retain good convergence characteristics. 
These stability considerations further emphasize the need to allow for a 
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reformulation of? the deterministic problem at a point sufficiently close to 
the design palirfc to obtain a good representation of the limit surface 
within the stdrtiHlty bounds imposed by the algorithm. 

2.LJ1.8 Inelastic Algorithm Development 

The iterative perturbation analysis algorithms developed 
for the linear eilastostatic problem Involve no assumptions on the linearity of 
the problem, and! conly minor coding and data management modifications should be 
necessary to extoend these algorithms to situations Involving mild 
nonlinearity- THie perturbation database must be extended to Include a record 
of the nodal strain histories, and the finite element code must be modified to 
carry along tnr {parallel the Incremental solution data for all perturbed 
problems. Perturbations must be allowed on additional types of variables, 
such as the material's elastoplastlc constants. This will Involve extensions 
to both the NESSUS/PRE and NESSUS/FEM modules. 

The extension of the perturbation algorithms In NESSUS/FEM to 
Inelastic problems raises important issues, which will affect the development 
of nonlinear algorithms for the remaining years of PFEM development. Version 
1.1 of NESSUS/FBM provides solution algorithms for deterministic linear 
problems using either a displacement-based or mixed Iterative finite element 
formulation. /KRll development of perturbation analysis algorithms to date has 
been based arr tfne displacement formulation. Implementation of the 
displacement metJhod for inelastic analysis will require changes to the 
internal data storage in the code, in order to retain the element strain 
history record aft the Integration points. All data Input and reporting of 
strains and stresses as perceived by the user can still be performed on a 
nodal basis. im alternate approach suggested by the NASA contract monitor 
Involves the edition of the mixed finite element formulation for 
probabilistic analysis, In order to maintain the node-oriented Internal data 
storage curreatTsy implemented in the code. This approach would lend Itself to 
a somewhat more elegant Implementation of some algorithms, but also involves 
substantial risfc associated with the adoption of less mature finite element 
technology. Before adopting the mixed approach for probabilistic finite 
element analysts.., the existing perturbation algorithms must be exercised and 
tested with simple elastostatlcs problems using the mixed formulation. This 
step is needed too ensure that no unexpected problems arise from the use of the 
current solution strategy with the mixed finite element formulation. If the 
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results of this experiment are positive, then It will be reasonable to 
consider proceeding with the mixed method for the development of Inelastic 
solution algorithms In NESSUS/FEH. 

Further development of Inelastic algorithms for probabilistic 
finite element analysis Is awaiting the outcome of the decision on the finite 
element formulation to be pursued. 

2.1.2 NESSUS/FPI Development 

2. 1.2.1 Eigenvalue Models for Non-normal Distributions 

Eigenvalue models for normally distributed, correlated 
variables (Reference [1]) have been used In the NESSUS code to solve problems 
Involving random fields such as pressure or temperature fields. The 
NESSUS/PRE module Is designed to generate uncorrelated variables based on 
the covariance matrix of the dependent variables. This Is described above 
In Section 2. 1.1. 3. The reasons for using the eigenvalue models are: 

1. Uncorrelated variables allow fast probability estimation 
using the NESSUS/FPI module. 

2. The number of the significant uncorrelated variables Is 
always less than the number of the correlated variables, and 
therefore the eigenvalue model Is able to reduce the dimension 
of the random variables entering the perturbation analysis. 

To extend the above method to problems Involving non- 
normally distributed, correlated variables, a model has also been formulated 
12 J. However, the new model Is still under Investigation and Is not yet 
Included In the NESSUS code system. A summary of the method Is given In this 
section. An example Involving a highly non-normally distributed variable 
Is given to test the model. The results suggest that. In order to obtain 
accurate transformed correlation coefficients, higher order terms In the 
Taylor's series expansion that Is used must be retained. Therefore, a more 
accurate formula relative to the one derived In [2] has been derived and is 
reviewed In this section. 

The eigenvalue model for the non-normally distributed, 
correlated variables requires two extra steps: the transformation of all 
variables to the normal distribution space; and, the derivation of the 
correlation coefficients of the transformed normally distributed variables. 

The normalization process Is defined in (9), 
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* W i * 1 » n (9) 

where F^(*) Is the original marginal distribution of the random variable 

x 1» where ♦() Is the normal CDF, and where Uj Is a standard normal 
variate. Note that (9) defines a one-to-one mapping; therefore, Xj may be 
formulated using the Inverse transformation: 


X 1 * F X i X (M u i)) (10) 

The Inverse CDF's, l.e., F^^ (•) are available In closed form for such 

distributions as the Welbull and Type I extreme value distributions. Using 
(10), the performance function becomes a function of u^. 

We next consider the computation of the correlation 
coefficients of the transformed variables. Consider two correlated random 
variables, denoted as Xj and X 2 . The correlation coefficients p„ w can be 
computed as 12 



E[X 1 X 2 1 - ECXjjEfXg] 
° 1°2 


( 11 ) 


Define the transformation from X^ to u^ as 


X 1 - T^(u^) 1=1,2 


( 12 ) 


and define 


H(u lt u 2 ) * XjX 2 = T 1 (u 1 )T 2 (u 2 ) (13) 

Eq. 11 may be expressed as 

°X 1 X 2 a l°2 * E -IH1 - EfTiUlV (14) 

Using a series expansion method, it can be shown that 

° X 1 X 2 V 2 * oIH ll + I ( h L 3 H 31> 1 + ?» 2h 22 + 


(15) 
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where higher order terms Is denoted (HOT) and where 


'Ij 


3ujau| 


d^T. 


u*0 


djl 2 
• IT 


duj du^ 


u=0 


(16) 


and where p Is the correlation coefficient of the transformed variables 
uj and U£. 

An approximation formula which Includes H^j terms up to 
1+J»8 has been derived earlier [21, but was later found to be Insufficiently 
accurate for very non-normal ly distributed functions, such as the uniform 
distribution. To Improve the accuracy of p, a more complete formula with 
up to 12 terms In the series has been derived and Is given In Table 2.2. 

In [2], a procedure was formulated to compute H^j using a 
numerical method. The procedure was demonstrated by using two examples 
Involving lognormal and normal variables. It was found that the series 
converges rapidly. However, It was not clear how the series would converge If 
the random variables were strongly non-normally distributed. In the following 
example a problem Involving a uniformly distributed variable Is tested to 
provide Information about the rate of convergence of the outlined 
procedure. This experience Is being used to guide the Implementation of 
the procedure Into the NESSUS code. 

Consider a case where one (say Xj) of the two random 
variables Is normally distributed; then 


dX 

du 


1 

1 


°1 


(17) 


d n X, 


du. 


= 0 for n > 1 


(18) 


for Uj = 0. 


Table 2.2 
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The approximating series for computing the transformed 
normal distribution correlation coefficient, p, may be derived as (up to the 
twelfth-order term): 


' 1 A 2 ‘ * k lHl1 + ? H 13 + 5 H 15 + k H 17 + 


p X,X o 0 2 


354 H 19 + 3P0 H l,ll 1 


(19) 


dx^ j d 3 x 2 
° du 2 ** 


u^ * 0 


where p Y Y Is the original correlation coefficient. 
a 1 a 2 

Assume that X 2 Is a unlformally distributed variable with 
a density function defined as 


f(X 2 ) = 1 
= 0 


0 < X 2 < 1 
otherwise 


( 20 ) 


Using the normalization scheme given In [1], the relationship between X 2 
and the standardized normal distribution function variable u 2 Is 

X 2 - #(u 2 ) (21) 

where *(u 2 ) Is the standard normal CDF. Eq. 21 Is plotted In Fig. 2.10. 

Note that a scale factor has been applied to X 2 such that a linear 
relationship with a slope of one in Fig. 2.10 represents a standard normal 
distribution. Therefore, the uniform distribution, according to Fig. 2.10 
behaves In a significantly non-normal ly distributed fashion for |u 2 | > 1. 

Using the numerical algorithm from [1J, thirteen sets of 
data are obtained as follows in Table 2.3: 



3 



Fig. MO Transforratlon from a Uniform Distribution to a Standardized Normal 



Table 2.3 


Results of Equivalent Normalization of Uniform Distribution 


Set 


u 2 


The next 


step 


1 

-6 

.990E-09 

2 

-5 

.287E-06 

3 

-4 

.316E-04 

4 

-3 

.00135 

5 

-2 

.00227 

6 

-1 

.158 

7 

-0 

.5 

8 

1 

.841 

9 

2 

.977 

10 

3 

.9986 

11 

4 

.999968 

12 

5 

.999999713 

13 

6 

.999999999 

Is to construct 

a twelfth-order polynomial 

denoted as 


X 


2 


12 

r 

n=0 


A u" 
n u 2 


( 22 ) 


The required derivatives for computing p are 




n! 


(23) 


where n = 1, 3, 5, 7, 9, and 11. 

By solving thirteen simultaneous linear equations, the 
coefficients A n can be found. Using Eq. 23 and Eq. 19 the approximation 
solution Is 


p y Y = 1.36225p[l. - 0.4380 + 0.2235 
* 1*2 


- 0.0871 + 0.0214 - 0.0024] 


( 24 ) 
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The final results are as follows 


p (Sixth-order) *0.9345 p Y Y 

* 1*2 

p (Eighth-order) * 1.051 p Y Y 

12 

(25) 

p (Tenth-order) * 1.020 p Y Y 

* 1*2 

p (Twelfth-order) * 1.023 p Y Y 

* 1*2 

The exact solution for this particular case Is available {31 and Is 

p (Exact) * 1.023 p^ x (26) 

Eqs. 25 and 26 show that the Taylor's series from (24) converges quite 
slowly. It needs to be pointed out that this example Is considered to be 
an extreme case to test the robustness of the algorithm. 

2. 1.2.2 FPI Validation Studies 

The original FPI (Fast Probability Integration) code 
using an algorithm developed by Wu [4J was modified to become the NESSUS/FPI 
code. In [4], the performance of the algorithm is assessed by six examples; 
some examples are considered the worst possible cases. The results Indicate 
that the algorithm Is able to provide accurate or reasonably good point 
probability estimates. In all cases, the results are significantly better 
than a widely-used FPI method: the first-order reliability analysis [5]. 

Chang [6] has Investigated the performance (accuracy and 
efficiency) of the FPI algorithm for computing structural reliability. 
Thirteen examples have been used to test the FPI accuracy. Many of the 
examples had nonlinear performance functions with non-normal random 
variables. The maximum number of random variables in the examples Is 
twenty. The results Indicate that the FPI algorithm provide good probability 
estimates. The errors In the point probability estimates are less than or 
near 5 % in twelve examples which are typical of mechanical design problems. 
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The only example which results In large error is the same 
problem Investigated In Reference (4). The performance function of the 
example Is a linear combination of the Identical and independent random 
variables. Each variable Is chi-square distributed with one degree of 
freedom. The distribution Is apparently highly non-normal. Its density 
function has a shape similar to an exponential density function (l.e., Aexp(- 
Ax) where A Is a positive constant) which has only one tall. It seems obvious 
that this distribution can not be fitted well by a symnetrlc bell-shaped 
normal curve. This test example Indicates that the accuracy of the current 
NESSUS/FPI code Is limited by the normality of the random variables. 

However, In probabilistic structural analysis, non-normal engineering 
variables are commonly modeled using the standard distributions such as the 
lognormal and the extreme value distributions. Using the FPI algorithm, 
these distributions can usually be fitted very well (In the least-squares 
sense) by the three-parameter normal distributions. 

The FPI code has also been compared with a code based on 
the second order reliability methods [7J. Three examples taken from Chang's 
report were tested. The comparison of computed probability estimates suggest 
that there are no significant differences In accuracy. The computational 
efficiencies were also compared by assuming that the computational 
efficiency for the first-order reliability analysis should be approximately 
equal using the two codes. 

The comparison of the computer time seems to confirm 
that, at least for linear performance functions, NESSUS/FPI Is faster than the 
second order reliability methods, especially for a large number of variables 
(in the test examples, the maximum number of random variables, N, Is 20). The 
reason Is believed to be that the second order methods needs to compute all 
the second order derivatives of the performance function In the transformed 
"standardized normal (u) space", whereas the NESSUS/FPI algorithm considers 
only part of the second order derivatives of the performance function in 
the X space. The advantage of using the NESSUS/FPI algorithm Is 
significant, since the computational effort required by the NESSUS/FPI Is 
of order N, while that required by the second order methods are of order 
N^. Moreover, since It Is very Inefficient to establish a "complete" 
quadratic response function In a typical NESSUS analysis, it seems more 
likely that the established response function will be either linear or 
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Incomplete quadratic. In such cases, the NESSUS/FPI algorithm Is 
particularly efficient because It does not require the computation of the 
second order derivative. 

The NESSUS/FPI algorithm has also been used to 
demonstrate how to compute the probability of Instability of a dynamic system 
[8J. The system Is represented by an n-th order linear differential 
equation. By assuming a solution of the form exp(st), a characteristic 
polynomial equation Is obtained where the coefficient are random functions of 
the random variables. A root s with a positive real part means that exp(st) 
becomes unbounded and the system Is unstable. A procedure based on the FPI 
algorithm Is developed and demonstrated using an example Involves a six 
degree polynomial with two random variables. For comparison purposes, a 
Honte Carlo solution Is obtained. The result shows that FPI is accurate 
and Is far more efficient than the simulation method. 

Other NESSUS/FPI validation exercises Include the 
solution of the NESSUS validation test cases 1 and 2 In which good agreement 
between FPI and Monte Carlo are obtained. 

A general conclusion drawn from the results of the 
numerous examples is that the NESSUS/FPI is consistently able to provide 
accurate results so long as the expansion point Is the most probable point. 
When the most probable point can be located (by Iteration), good results can 
usually be expected even with linear approximation of the performance 
function. 

2. 1.2.3 FPI Accuracy/ Improvement Studies 

A number of studies on the FPI algorithm were conducted 
at the University of Arizona. The studies focused on approximation functions 
within NESSUS/FPI which have been suspected of producing errors In the 
resulting probability estimates. Modifications to NESSUS/FPI have been made 
to Improve the performance of the code and Include: 

1. A new gamma function has been Introduced. This function 
representation has about eight significant figures for accuracy 
and Is a significant Improvement over the polynomial 
approximation previously used. 

2. Changes were made In calculations of the extreme value 
distribution (EVO) parameters to provide ten-place accuracy. 
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3. The polynomial approximation to the Inverse normal CDF has 
been replaced by the secant method with a significant 
Improvement In accuracy. 

4. All distribution parameters are now computed at the 
beginning of the program Instead of within the subroutines. 

5. The secant method Is used to compute the Uelbull shape 
parameter. Using this method and the new gamma function should 
Improve accuracy of both Welbull parameters. 

Details of these changes are reported In Appendix C. 

Comparison of the old and new code for HESSUS/FPI showed 
small changes In the results, generally less than 5 % for all thirteen test 
examples. However, the new code accuracy has been achieved with no 
significant loss of efficiency and Is, therefore, being Incorporated In the 
next release of NESSUS/FPI. 

In addition to the above numerical Improvements, there 
are three new distributions which have been added to MESSUS/FPI. The new 
distributions are: 

1. The Frechet distribution Type 2 asymptotic distribution of 
extreme values from an Initial lognormal distribution and, In 
general, an Initial distribution having a polynomial tall In the 
direction of the extreme. 

2. Truncated Welbull distribution. 

3. Truncated normal distribution. 

The truncated distributions are included for modeling distributions of 
material axes for the turbine blade verification problem. Other 
distributions already in the code are the normal, lognormal, Welbull, Type 
1 extreme value, maximum entropy, chi-square, and NESSUS. The NESSUS 
distribution Is a polynomial of a normally distributed variable. 

2. 1.2.4 Confidence Band Estimation 

The basic goal of confidence band estimation. In the 
context of the NESSUS analysis. Is to quantify the confidence on the 
accuracies of the probability estimates for the response functions. The basic 
assumption for the methods Is that the response functions are derived from the 
NESSUS perturbation data base. The approach Is to treat the distribution 
parameters of the input random variables as random variables, and then 
create the COF of the response function CDF. This strategy Is the essence 
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of the Beyslan approach to parameter estimation. Four approximation 
methods are Identified for estimating the confidence (or error) band of the 
cumulative distribution function of the response function. The most 
suitable one Is Identified and Included In the new version of the 
NESSUS/FPI code. 

Let Z denote the response variable, and Z(X) denote the 
response function where X Is a vector of the Input Independent variables. In 
the NESSUS/FPI, Z(X) is a polynomial function: 

N N 

Z - Z(X) * a 0 + r a 1 x 1 + z b.X? (27) 

where N Is the number of Independent random variables. In general, several 
polynomial equations may be required to ensure sufficient accuracy of the 
function over a wide range of Z. Ideally, one polynomial should be 
established for a selected Z value. 

The basic assumption In the response function for the 
confidence Interval estimation Is that for a given Z, the best estimate Z(X) 
(derived using the best estimates statistics of X) remains valid within the 
confidence band. In general, Z(X) Is different 1 for different distribution 
parameters set because the most probable point, which Is used to define 
Z(X) , Is a function of the distributions. However, the assumption Is valid 
when Z(X) Is actually a first or second degree polynomial. For highly 
nonlinear Z(X) function, the assumption is a reasonable one so long as the 
variabilities of the significant random variables are not very large. 

There are two basic types of uncertainties in a 
NESSUS/FPI-generated response function: (1) physical uncertainty and (2) model 
uncertainty. Physical uncertainty Is the uncertainty associated with physical 
phenomena which are Inherently random. In the NESSUS analysis, this 
uncertainty Is accounted for by treating the Input variables as random 
variables or random fields. Model uncertainty Includes parameter uncertainty, 
uncertainty In the statistical distribution model, response function model 
error, etc. The approach adopted In this study concentrates on the 
variabilities of the Input variables. 

Assume that X is a normally distributed random variable 
with mean y and standard o deviation. Given a sample, n, the sample mean, Y, 
Is a normal variable with mean and standard deviation of 



requires major code modification and developement effort. Method 3 Is 
accurate for large samples, but the full simulation Is extremely 
Inefficient. Method 4 Is accurate and Is consistent with the current 
NESSUS/FPI approach, l.e., the CDF of Z(X) Is computed by using the FPI 
method for a given set of statistical parameters. In terms of the 
computational efficiency. Method 4 Is Inefficient relative to Methods 1 and 
2 but Is much faster than Method 3. Overall, Method 4 was considered 
accurate with satisfactory efficiency, therefore. It was selected and 
has been Incorporated In the NESSUS/FPI code. 

An example has been taken from that proposed In Appendix D as a means 
to test the confidence band estimation algorithm In the NESSUS/FPI code. 

The response function Z Is a function of two normal variables X and Y, 

Z = X - Y 

The statistics are. 

For X For Y 

n = 20 n = 20 

Y = 10 Y - 5 

S)( * 2 Sy * 1 

The sample means and the sample standard deviations are defined as the best 
estimates. Using Eqns. 28 and 29, the COVs for the means of both X and Y 
are 0.0447; the COVs for the standard deviations of both X and Y are 0.162. 

By entering these parameters into the NESSUS/FPI, the 90# and 95X 
confidence bands of the CDF of Z were generated. The result Is shown In 
Fig. 2.11. The Monte Carlo sample size Is 5,000. 

2. 1.2. 5 Monte Carlo Methods 

Monte Carlo simulation has been usually considered to be 
last resort for solving a major simulation problem because of Its high cost 
for accurate results, especially In the tails of the distributions. However, 
recent developments of new and efficient algorithms have made Monte Carlo 
more attractive. 

A study of several Monte Carlo simulation algorithms has 
been conducted at the University of Arizona for the PSAM project. Two 



Fig. 2.11 Confidence Interval Estimation 
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^ « u s » a// n (28) 

2 2 

The variable (n-1) s /o has a chi-square distribution with n-1 degree of 
freedom. If this chi-square distribution Is approximated by a lognormal, then 
the distribution of s will be a lognormal. The statistic of s can be 
approximated as 

u s * s o $ = a// [ 2(n-l) ] (29) 

For the NESSUS/FPI confidence band estimation, we assume 
that each X^(1=1,N) Is characterized by Its mean and standard deviation. We 

further assume that the statistical distribution Is normal for the mean, and 

Is lognormal for the sample standard deviation. These assumptions about the 
statistical distributions of the parameters are exact only when X Is normal. 
The actual distributions usually do not follow available standard 
distributions and the CDF's cannot be defined In closed forms. 

The required Input data for the confidence band 
estimation are the statistics (the means and the COVs (coefficient of 
variation « standard deviation/mean)) for the means and the standard 
deviations of all X^'s. Note that the input statistics may be estimated by 
using Eqs. 28 and 29 where the actual statistics may be replaced by sample 
statistics. However, the statistics may also be estimated using other 
statistical methods or engineering judgement. This Input format Is more 
flexible since It does not require that the sample sizes be defined. However, 
the Input statistics must be prepared before the estimation process. 

Four methods are considered: 

1. First Order Error Bounds 

Assume that Z(X) Is linear and each X 1 Is normally distributed. For 
each Z(X) , there is a best-estimate most probable point of X. The best 
estimate CDF of Xj, denoted as F^, Is determined using the most probable 
point for each X^. At F^, X^ can be written In terms of the mean and 
standard deviation by Inverting the CDF, 

X i 3 F i J + ^ 


(30) 
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where and are random variables. Upon substitution of each Into Eq. 
30, Z(X) can be expressed In terms of Z(y,o). , 

By further assuming that Z(X) has a normal or a lognormal 
distribution, closed form solutions for the confidence bounds were derived 
by Wlrschlng and are Included In Appendix D. 

2. FPI algorithm 

Assume that Z(X) Is linear. The CDF of Z(X) = z, In terms of the 
standard normal variate, u, can be formulated as 


a Q + ra i y 1 

7 ~ 9 — 

/ (za^) 


(31) 


where and are the equivalent normal parameters of X^ based on the 
Rackwltz-Flessler algorithm [9]. Eq. 5 Is a safety Index formulation 
based on the first order reliability method. 

Note that the equivalent normal parameters are functions of the CDF, 
F(X), the POF (probability density function), f(X), and the most probable 
point X.j . Let 


9 i = (Vj.Oj) 


(32) 


Because F(x) and f(x) are functions of *, therefore, u can be expressed as 


u = function (e) (33) 

The limit state or performance function can be formulated as 

9(?) = u - u Q (34) 

The following Is a proposed FPI Iteration algorithm for estimating the CDF of 
u for a selected response function value Z * z: 

1. Select a u Q . 

2. Guess the design point of the distribution parameters, e. 

3. Compute the equivalent normal parameters of the random 
variables, e. 


45 


4. Define the distributions of X using the most probable point 
of 8. 

5. Guess the most probable point of the basic variable X. 

6. Compute equivalent normal parameters for non-normal X. 

7. Compute the most probable point of X and the CDF of Z(X) * z 

8. Go to step 3; repeat until the most probable point of X or 
the CDF of Z(X) Is stabilized. 

9. Compute the most probable point of e and the CDF of g (e) - 0 

10. Go to step 2; repeat until the most probable point of e or 

the CDF of g(e) = 0 Is stabilized. 

Note that the above procedure requires nested iteration loops. Step 3 to 
step 8 constitute the Inner FPI loop for a selected e set. Steps 2, 9 and 
10 constitute the outer loop for finding the most likely e set. 

3. "Full 11 Monte Carlo Simulation 

This method is conceptually more straight-forward. It requires the 
following steps: 

1. Generate samples of e sets, e., j * 1, J 

2. For each e, generate a set of X k , k = 1, K 

3. Compute, using X., the response function value, Z^, k = 1,K 

4. For each 6j, compute the CDF of Z(X)=z, denoted as (CDF) j , 
j = 1, J 

5. Using samples of (CDF), construct COF of u. 

This last procedure is expected to be extremely time-consuming because It 
requires the generation of "J times K" samples of Z(X) values. 

4. Mixed Approach - Combination of Monte Carlo and FPI 
This approach combines the above methods (2) and (3). The difference 
between this approach and the previous approach (Method 3) Is that after a 
set of Xj Is generated, the FPI routine is used to compute each (CDF) ^ . 

The methods can now be compared. Method 1 captures the essence of 
statistical uncertainty and is the most efficient. However, the accuracy 
of Method 1 Is limited by the distributional assumptions. Further 
improvement is needed for this fast algorithm. Method 2 has the potential 
to be both fast and accurate, however, it is the most complicated and 
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computer programs based on the conventional Monte Carlo algorithm and one 
based on variance reduction using antithetic variables were written at the 
University of Arizona. Other efficient Monte Carlo schemes are still being 
evaluated. The work to date Is summarized In Appendix E. 

2. 1.2.6 Integration with NESSUS/FEM 

This section summarizes the study of an FPI Iteration 
procedure which was Intended to be used to Integrate the FPI and the FEM 
modules. The procedure was used successfully to solve several selected 
problems. At the end period of this study, however, a new and potentially 
more efficient method was formulated which seems to be most suitable for 
constructing the CDF of a response function. The newly-developed method and 
the Iteration procedure are summarized In the next section (2. 1.2. 7). The 
procedure described in the present section Is useful for computing a point 
CDF. For creating the entire CDF, the present procedure may ultimately be 
replaced by the new procedure. However, the new procedure is based on the 
present study, and many key concepts discussed In this section remain valid 
for the new procedure. 

The integration of the NESSUS/FPI and the HESSUS/FEM Is 
based on the concept of successive 11 near /quadrat 1c approximation algorithm 
Which was identified In the first year of this project [5J. The goal Is to 
expand or perturb the performance function about the most probable point. 

Mote that In the field of structural reliability analysis, where the goal Is 
to find the probability of failure estimate, the most probable point Is called 
the "design point". The algorithm which Is summarized In the following has 
been used to test several examples with success. 

The Iterative algorithm has been established as follows: 

1. Identify critical dependent variable (stress, frequency,...) 

2. Select values for dependent variable, (e.g., mean, mean + 

103» of mean) 

3. Using the NESSUS/FEM module, compute the perturbation 
solutions about an Initial guessing most probable point. 

Initially this can be chosen as the mean values. However, a 
good Initial guess of the most probable point will accelerate 
the Iteration procedure. 

4. Establish linear/quadratic response surface from the 
perturbation solutions using the least-squares method. 
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5. Compute the most probable point, using the NESSUS/FPI, for 
the selected value of dependent variable. 

6. Compute the CDF and evaluate accuracy (based on the 
successive CDF values or the most probable point values). 

7. Use the new most probable point and go to step 3, until the 
solution converges. 

Experience with this algorithm has Indicated that the solution can usually 
be found In about three Iterations. 

An example Is now presented to Illustrate the above 
Iteration procedure for Integrating FEM with FPI. The example Is a simplified 
version of the NESSUS validation test case 2 from the first year annual 
report. The finite element model Is Illustrated In Fig. 2.12. There were 
Initially ten random variables: five correlated loadings, width, length, 
thickness, base spring and modulus of elasticity. By assuming that the width 
Is deterministic, the "exact" root stress becomes: 

S * LP/t Z (35) 

In which P Is a load random variable; L Is the length and t Is the 
thickness. The mean value of S Is approximately 3500 psl. 

In order to illustrate more clearly the Iteration 
procedure, It is assumed further that t Is a deterministic variable and L an P 
are normally distributed. Note that none of the above assumptions Is 
required for the NESSUS solution. 

Define the "reduced variables" of L and P as 

U 1 * *-mean^ L std. dev. 

(36) 

u 2 * (P - p mean^ p std. dev. 

Using Eq. 36, L and P can be expressed as functions of Uj and 
respectively. Substituting L and P Into (35), one can plot the contours of 
constant stress (Iso-stress)in a two dimensional u space as shown In Fig. 

2.13. The reason for using the u coordinate system Is that the joint 
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Fig. 2.12 The Test Example 
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Fig. 2.13 First Iteration Used to Obtain the Most Probable Point 
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probability density function Is rotatlonally symmetric. The most probable 
point is, therefore, easily Identified as the point on a Iso-stress curve 
which Is nearest to the origin. 

The problem Is defined as follows: Find the most probable point (and 
the COF) for S 2 * 4500 psl by starting at Sj » 3500 psl. Fig. 2.13 shows 
the result of the first Iteration. Initially the linear approximation of S 
Is based on the mean values of P and L which corresponds to the origin. A 
■predicted" Iso-stress curve (S 2 * 4500 psl) can be defined using the 
mean-derived linear equation. The predicted S 2 curve, which Is parallel to 
the approximated Sj curve, deviates from the exact S 2 curve because S Is 
actually a nonlinear function of P and L. However, this first Iteration 
leads one to a region close to the exact most probable point. Using the 
predicted most probable point as a new expansion point, a second Iteration 
results In an accurate prediction of the most probable point as shown In 
Fig. 2.14. No more Iteration Is required. 

For S 2 > 4500 psl, the volume under the joint probability density 
function surface Is concentrated near the most probable point. The first 
order reliability analysis gives the following result: 

P(S > 4500 psl) = ♦(- 6) (37) 

where 6 Is the minimum distance defined by the most probable point. 

The above procedure can be applied to several values of S In order to 
eatabllsh the entire CDF. In the following, the procedure will be applied 
to Integrate the FPI and the FEM. The test problem Is the NESSUS 
validation test case two of which the width Is assumed to be deterministic. 
The results of the first Iterations at three stress levels (2600, 3500 and 
5400 psl) are shown In Fig. 2.15. 

The purpose of Fig. 2.15 Is to show the algorithm for Integrating the 
NESSUS/FEM and the NESSUS/FPI. The finite element model consists of twenty 
plate elements (NESSUS element 75). The difference between the analytical 
and the NESSUS solutions is about 3%. In order to show the effect of 
successive linear approximation, a "calibrated exact" CDF Is used to match 
the mean solutions. 

The first perturbation was taken about the mean values of the 
independent variables. Two more FEM perturbations were taken about S = 


’New 1 predicted S 
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Fig. 2.14 Second Iteration Used to Obtain the Most Probable Point 
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Fig. 2.15 Example - NESSUS/FEM and NESSUS/FPI Integration {Stress) 
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2600 and S * 3500 using the predicted most probable points. It Is shown 
that the CDF values are accurate "locally" around the approximation 
points. 

Figure 2.16 shows the result of an analysis for the tip displacement of 
the same cantilever beam. The goal was to compute the CDF at 1.2 Inches. 

The result of the first Iteration Indicates a significant Improvement at 
the region around 1.2 Inches. 

It should be noted the above results were obtained using "small" 
perturbations (0.05 or 0.1 standard deviation for the Independent random 
variables). The reason was to estimate the first order sensitivities more 
efficiently. It Is noted also that the update of the most probable points 
In the NESSUS/FEM Input data deck were done manually. The updated 
"correlated" nodal loads were being computed using the most probable point 
values of the "uncorrelated" loads (which means that the eigenvectors 
generated using the NESSUS/PRE module must be used to update the NESSUS/FEM 
data). This computational procedure needs to be considered carefully In 
designing the user-friendly expert system - the NESSUS/EXPERT module. 

2. 1.2.7 A New CDF Estimation Method 

This section summarizes a new CDF estimation 
method. This method, If proved to be more effective for estimating the CDF of 
the response function, will replace the one described In the previous section 
(2. 1.2.6). However, since the new method was developed In the last period of 
the second year PSAM efforts, further detailed study of the method Is required 
to validate the method. A preliminary discussion on the method Is given In 
Appendix F where the formulation of the method and a simple example are 
Included. By using a procedure which corrects the error of the response 
function at the most probable point, It Is shown that the new procedure has 
the potential to significantly Improve the NESSUS solution efficiency by 
reducing the requirements on the perturbation solutions. 

The procedure based on the new method for integrating the 
NESSUS/FEM and the NESSUS/FPI modules Is as follows: 

1. Construct first (can be extended to second) order 
approximation of the response function Z(X) about the mean 
values. NESSUS/FEM module Is used! to generate response 
function sensitivities or perturbation solutions. 

2. The reponse function Is established using the least-squares 
routine In the NESSUS/FPI. 
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3. Using the response function of step 2, a COF can be 
constructed. (This COF Is, In general, not sufficiently 
accurate at the tall regions of the distribution.) 

4. Select a CDF value from the result of step 3; find the 
corresponding "predicted" response value, Zp. 

5. At Z * Zp, compute the most probable values of X, Xp, using 
the NESSUS/FPI module. 

6. Using the NESSUS/FEM, compute the "exact" response function 
value, Z e , at the most probable point, X p . Z e Is the 

"corrected" value for the selected COF value defined In step 4. 

7. Compare Z e and Z p . If the difference Is small (say, less 

than 20 X) go to step 3 and select another probability level. 

If the difference Is large, go to step 1 and replace the mean 
values of X by the X p values. 

The significant difference between the present procedure and the one 
presented in the previous section is that the present procedure fixes a CDF 
value and looks for the accurate corresponding response function value, 
whereas In the previous procedure, a response function value Is fixed and 
the CDF value Is found using an iteration procedure. Thus, the present 
procedure relies more on the additional deterministic solutions while the 
previous procedure relies heavily on the additional sensitivity analyses. 

Since the sensitivity analyses require more computational efforts than the 
deterministic analyses, it seems resonable to expect that the new procedure 
will be more efficient. 

2. 1.2.8 NESSUS/FPI Code Validation Studies 

A test plan for validating the first year 
probabilistic finite element code was Included In the First Year Annual Report 
(Vol. Ill, Section 4). It consisted of nine validation problems which were 
designed to test a variety of capabilities of the NESSUS code. The exact 
solutions, In terms of the probability distributions or the probability of 
exceedance, have been obtained for the first five validation problems. The 
results which are summarized in the following are presented In a format 
compatible with the NESSUS/FPI output. "Exact" solutions are obtained using 
the Monte Carlo simulation if no closed form solution Is available. These 
exact solutions are to be compared with the NESSUS solution to validate the 
code as well as the solution algorithm (i.e., FPI iteration algorithm). 
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The exact solution for the validation problem 1 was 
Included In the First Year Annual Report. The 'problem addressed Is a 
cantilever beam subjected to static loadings, P^I-1,5) (see Fig. 2.17). The 
loadings are correlated random variables. Other random variables Include 
Young's modulus, length, thickness, width, base spring and yield strength. 

The expected output of the NESSUS solution include the COFs of the maximum 
stress and the tip displacement, and the probability that the stress will 
exceed the yield strength. The type of finite element used In NESSUS Is 
beam element 98. 

Problem 2 Is similar to problem 1 except that the plate 
element is used and the thickness of the beam Is reduced. Because of the 
reduced thickness, the nodal loads were changed from 20 lbs to 0.1 lbs. Figs. 
2.18 and 2.19 summarize the results of the COF of the maximum stress, the CDF 
of the tip displacement and the probability that the stress exceeds the yield 
strength. 

The goal of the validation problem 3 Is to validate the 
NESSUS eigenvalue solution algorithms. The cantilever beam defined In problem 
1 Is used again. The response functions tested are the first three bending 
frequencies In each lateral direction. The analytical solutions for the 
frequencies In the X and Z directions modes were used to derive the exact 
CDFs. Figure 2.20 and Table 2.4 summarize the results for the X direction; 
Fig. 2.21 and Table 2.5 summarize the result for the Z direction. 

Validation problems 4 and 5 addressed a rotating beam as 
Illustrated In Fig. 2.22. The random variables are: mass density, length. 
Young's modulus, thickness and width. Problem 4 tests the beam element, and 
problem 5 tests the plate element. The response function tested are the 
tip axial displacement and the first bending frequency In the Z direction. 

The analytical solutions are the same for both problems. In the original 
test plan, the beam was fixed at the rotation center. To represent a 
turbine blade configuration more closely, the Inner radius (measured from 
the center of rotation to the "fixed" end of the beam) was defined to be 
4.237 inches. Analytical solutions were revised and used to generate exact 
solutions using Monte Carlo simulation. Figures 2.23 and 2.24 summarize the 
results for the tip displacement and the fundamental bending frequency. 

The NESSUS code validation Is still In progress, and MARC 
will run the NESSUS/FPI code and compare results with these "exact" 
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Fig. 2.17 Model Definition of a Cantilever Beam 
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Fig. 2.19 Validation Test Case 2 

(Exact COF of Max. Stress) 
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Fig. 2.20 Validation Test Case 3 

(Exact CDF's of Frequencies (+/-X Dir.) 
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Table 2.4 


Validation Case 3 Cantilever Beam (Natural Frequency) 
X Dir. (Horizontal) 



1st Mode 
2nd Mode 
3rd Mode 

Median 

508.04 

3233 

8905.16 

Mean 

511.4438 

3254.661 

8964.824 

Std. 

59.28826 

377.2911 

1039,232 


w(rad/sec) 

U 1 

u 2 

w(rad/sec) u 2 

u 3 


300 

-4.56199 


2567.145 

-1.99717 


315 

-4.13945 


2695.502 

-1.57464 


330.75 

-3.71692 


2830.277 

-1.15211 


347.2875 

-3.29439 


2971.791 

-0.72958 


364.6518 

-2.87186 


3120.380 

-0.30705 


382.8844 

-2.44933 


3276.399 

0.115481 


402.0286 

-2.02680 


3440.219 

0.538013 


422.1301 

-1.60426 


3612.230 

0.960544 


443.2366 

-1.18173 


3792.842 

1.383076 


465.3984 

-0.75920 


3982.484 

1.805608 


488.6683 

-0.33667 


4181.608 

2.228140 


513.1018 

0.085857 


4390.689 

2.650671 


538.7568 

0.508389 


4610.223 

3.073203 


565.6947 

0.930921 


4840.734 

3.495735 


593.9794 

1.353453 


5082.771 

3.918267 

-4.85640 

623.6784 

1.775984 


5336.910 

4.340798 

-4.43387 

654.8623 

2.198516 


5603.755 

4.763330 

-4.01134 

687.6054 

2.621048 


5883.943 


-3.58881 

721.9857 

3.043580 


6178.140 


-3.16627 

758.0850 

3.466111 


6487.047 


-2.74374 

795.9893 

3.888643 


6811.400 


-2.32121 

835.7887 

4.311175 


7151.970 


-1.89868 

877.5782 

4.733707 


7509.568 


-1.47615 

1824.422 


-4.95489 

7885.047 


-1.05362 

1915.643 


-4.53236 

8279.299 


-0.63108 

2 011.425 


-4.10983 

8693.264 


-0.20855 

2111.996 


-3.68730 

9127.927 


0.213974 

2217.596 


-3.26477 

9584.324 


0.636505 

2328.476 


-2.84224 

10063.54 


1.059037 

2444.900 


-2.41970 

10566.71 


1.481569 




11095.05 


1.904101 




11649.80 


2.326632 




12232.29 


2.749164 




12843.91 


3.171696 




13486.10 


3.594228 




14160.41 


4.016759 




14868.43 


4.439291 




15611.85 


4.861823 
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Fig. 2.21 Validation Test Case 3 
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Table 2.5 

Validation Case 3 Cantilever Beam (Natural Frequency) 
Z Dir. (Horizontal) 


Median Mean Std 



1st Mode 
2nd Mode 
3rd Mode 

497.9 

3168.5 

8727.4 

501.2359 

3189.728 

8785.873 

58.10493 

369.7639 

1018.487 


w(rad/sec) 

U 1 

u 2 

w(red/sec) u 2 

u 3 


300 

-4.38739 


2830.277 

-0.97759 


315 

-3.96486 


2971.791 

-0.55506 


330.75 

-3.54233 


3120.380 

-0.13252 


347.2875 

-3.11979 


3276.399 

0.290003 


364.6518 

-2.69726 


3440.219 

0.712534 


382.8844 

-2.27473 


3612.230 

1.135066 


402.0286 

-1.85220 


3792.842 

1.557598 


422.1301 

-1.42967 


3982.484 

1.980130 


443.2366 

-1.00714 


4181.608 

2.402661 


465.3984 

-0.58460 


4390.689 

2.825193 


488.6683 

-0.16207 


4610.223 

3.247725 


513.1018 

0.260455 


4840.734 

3.670257 


538.7568 

0.682986 


5082.771 

4.092788 

-4.68178 

565.6947 

1.105518 


5336.910 

4.515320 

-4.25925 

593.9794 

1.528050 


5603.755 

4.937852 

-3.83672 

623.6784 

1.950582 


5883.943 


-3.41419 

654.8623 

2.373114 


6178.140 


-2.99166 

687.6054 

2.795645 


6487.047 


-2.56912 

712.9857 

3.218177 


6811.400 


-2.14659 

758.0850 

3.640709 


7151.970 


-1.72406 

795.9893 

4.063241 


7509.568 


-1.30153 

835.7887 

4.485772 


7885.047 


-0.87900 

877.5782 

4.908304 


8279.299 


-0.45647 

1824.422 


-4.78037 

8693.264 


-0.03393 

1915.643 


-4.35784 

9127.927 


0.388592 

2011.425 


-3.93531 

9584.324 


0.811124 

2111.996 


-3.51278 

10063.54 


1.233656 

2217.596 


-3.09025 

10566.71 


1.656187 

2328.476 


-2.66771 

11095.05 


2.078719 

2444.900 


-2.24518 

11649.80 


2.501251 

2567.145 


-1.82265 

12232.29 


2.923783 

2695.502 


-1.40012 

12843.91 


3.346314 




13486.10 


3.768846 




14160.41 


4.191378 




14868.43 


4.613910 




15611.85 


5.036441 
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Fig. 2.23 Validation Test Cases 4 & 5 
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Fig. 2.24 Validation Test Cases 4 & 5 
(Exact CDF of 1st Frequency) 
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solutions. To date these validation problems have been used to uncover 
several program bugs, to gain experiences for Incorporating user-friendly 
Interfaces, and to lead to new solution strategies. 

The validation of problem 1 has raised an Issue regarding 
the random variables data Input structure. The type of finite element used In 
this problem Is Timoshenko beam element (NESSUS element type 98). The random 
variables Include thickness, t, and width, w, among others. The first-year 
NESSUS code defines the beam section using the area. A, and the area moment 
of Inertias I x and I y . This format Is not proper because A, I x and I y are 
correlated depending on the shape of the beam sections. Conseqently, the 
Independent perturbations of w and t are Impossible. To correct the 
dependency problems requires that the NESSUS/FEN code use "basic variables" 
w and t as input data. This strategy can be applied to other problems 
Involving dependent variables. 

Pending implementation of t and w as random variables, 
problem 1, with w and t as deterministic values, was used to validate other 
capabilities of the code. Modal frequencies, stress and displacement 
solutions were obtained and compared well with the analytical solutions. The 
perturbation solutions were not obtained, however, pending the code 
modification of the input structure. 

For the validation problem 2, perturbation convergence 
Instability has been observed for the width, w. In order to obtain a complete 
perturbation data base and to accelerate the validation process, w was 
temporarily treated as a deterministic value. The validation study of this 
slightly modified problem 2 has resulted in the successful integration of the 
NESSUS modules (PRE.FEM and FPI), using a successive linear approximation 
algorithm (Section 2. 1.2.6). The study of the FPI Iteration procedure for 
this problem has also led to a new and potentially more efficient solution 
strategy for estimating CDF (Section 2. 1.2.7). 

A validation problem not Included In the first-year plan 
Is a simple model simulating a turbine blade. The goal Is to validate the 
capability of the code to treat the material axes as random variables. The 
model consists of four solid elements (NESSUS element type 75). The material 
has anisotropic property with one material axis modeled as a random variable. 
Perturbation results for the first two modal frequencies were obtained to 
estimate the CDFs using the NESSUS/FPI. Analytical solutions for this test 
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case Is not available. However, the validation results were judged 
reasonable based on the Information of the resulting means and standard 
deviations of the natural frequencies. A data Input limitation was 
Identified In that the material property (0) matrix defined In the FEM 
Input data Is deterministic. That Is, the material properties such as the 
Young's modulus and the Poisson's ratio cannot be defined as random 
variables. It appears that code modification Is necessary to solve the 
problem. — r 

MARC has now completed the perturbation analysis for the 
validation problems 1, 3 and 5. New CDF estimation procedure (Section 
2. 1.2.7) will be used to continue the validation of the NESSUS modules and the 
solution procedure. 

2.1.3 NESSUS/EXPERT Development 
2. 1.3.1 Approach 

The goal of the expert user Interface Is to provide a 
flexible, user-friendly interface to the NESSUS/FEM and NESSUS/FPI codes. 

This Interface will serve not only as an enhanced, on-line, automated user's 
manual for these codes, but It will also act as an expert aid In generating 
a data deck for a problem, especially the probabilistic Information needed 
to solve a problem using NESSUS. Emphasis has been placed on minimizing 
the detailed knowledge that a user must have of NESSUS, allowing him/her to 
provide the Information about a particular problem In as natural a way as 
possible and having the the expert user interface generate the actual data 
deck required. 

To this end, an expert system called NESSUS/EXPERT Is 
under development in parallel with the development of the NESSUS code Itself. 
The system will consist of two major components, the Interface to NESSUS/FEM 
and the Interface to NESSUS/FPI. The Interface to NESSUS/FEM is to contain 
essentially all of the knowledge about the use of NESSUS provided In the 
user's manual. It will also contain any clarification and other specifics 
about the use of the code known to those who developed the code and those who 
have tested it. It will also contain knowledge about generating probabilistic 
Information about the problem from general descriptions. The Interface to 
NESSUS/FPI will contain knowledge on how to analyze and interpret the results 
of a run, thus aiding the user in deciding what to do next. 
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Host such knowledge Is embodied In the form of rules-of- 
thumb that provide methods for calculating specific values needed by HESSUS 
given general Information about the problem, that provide hints about how best 
to use the system, and that Indicate what Is useful and Important In the 
output of a run. Thus, the problem fits. In a fairly straightforward 
manner, the production rule knowledge representation method. This Is 
convenient since most existing expert system building aids are rule-based 
and this Is the best understood method of the AI technologies. Thus, the 
approach Is to design and Implement two rule-based expert systems to act as 
an Intelligent front and back end to the NESSUS code. 

2. 1.3.2 LISP/0PS5 Environment 

The programming language selected for Initial 
development of NESSUS/EXPERT Is 0PS5. 0PS5 Is an expert system building 
software facility that allows a programmer to write production rules directly 
as code. The version of 0PS5 being used In NESSUS/EXPERT Is public domain and 
available free from Carnegle-Mellon University. It runs under Franz Lisp Unix 
on a OEC VAX. SwRI has recently ported this version of 0PS5 to DEC Common 
Lisp so that It now runs under DEC VMS and on the Sun Workstation under Sun 
Common Lisp. 

The entire NESSUS/EXPERT system will be coded Inti ally 
using 0PS5. The advantage of such a tool Is that It offers a much higher 
level of productivity for the programmer because the knowledge can, to some 
extent, be encoded directly Into 0PS5 code. It also produces a much more 
readable and maintainable computer program. Though there are many other more 
elaborate, and more expensive, methods and tools for creating expert 
systems, the production rule technology embodied In 0PS5 Is sufficient for 
this task. 

The major drawback of using a tool such as 0PS5 for this 
application Is its dependence on the Lisp environment. 0PS5 Is an interpreter 
coded In Lisp and, therefore, requires Lisp In order to run. Lisp does not 
currently provide an easy Interface to FORTRAN on the DEC VAX. Thus, In the 
case where a data deck Is produced for the pre-processor, the pre-processor 
cannot be Invoked directly from NESSUS/EXPERT. Instead, the user must 
leave NESSUS/EXPERT, run the FORTRAN-based pre-processor, and then return 
to NESSUS/EXPERT where the resulting file can be read In and the process of 
developing a data deck for NESSUS can continue. 
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A Lisp-based system Is being used because there are 
currently only about three expert system building tools available written In 
FORTRAN. Their functionality Is limited and the resulting code Is not all 
that readable because of the compromises made due to the FORTRAN language. A 
solution to this problem, as well as the requirement that all code for the 
PSAM project be delivered In FORTRAN, is to reimplement 0PS5 In FORTRAN. 

The dependence on the Lisp environment would be removed and the Interface 
to FORTRAN would be automatic. Another option would be to recode the 
entire NESSUS/EXPERT system In FORTRAN at the completion of this project. 

This Is not desirable because all of the flexibility and maintainability 
acquired by using 0PS5 will be lost In the translation. Therefore, for the 
moment NESSUS/EXPERT will remain In Lisp-based 0PS5. 

2. 1.3. 3 NESSUS/FEM Interface 

Development of NESSUS/EXPERT has begun with the 
creation of an expert system for Interfacing to NESSUS/FEM. Because the 
expert system developed must be an "expert" In how to use NESSUS/FEM, work has 
started by Incorporating the knowledge contained In the MHOST User's Manual. 
Examination of the MHOST User's Manual supplied by MARC Analysis Resesarch 
Corporation has revealed a list of various types 1 of knowledge that must be 
•used when creating a data file for NESSUS/FEM that will run correctly for a 
specific problem. These Include: 

1. The required information for all problems ( 1 .e. , number of 
elements, connectivity of the nodes, etc.) 

2. Interdependencies of options selected and data provided with 
other possible options and data (l.e., the number of elements 
provided under the model data must be less than or equal to the 
number provided under the parameter data, the *compos1te option 
under parameter data requires the ‘laminate option under the model 
data, etc.) 

3. Incompatible selections of options/data (l.e., ‘frontal solution 
option cannot be used with the ‘bandmatrlx option) 

4. Default options and values (l.e., ‘bandmatrlx is the system 
default option, upper bound to the number of beam element crossings 
defaults to 1, etc.) 

5. All available keywords and their "meanings" 

6. Format of the parameters and data expected for each keyword, both 
for acquiring the needed information and for setting up the data 
file properly 
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7. Helpful hints concerning Idiosyncrasies of the code (l.e., not 
recommended to use the option ‘stress, not recommended to use the 
option ‘displacement method for Inelastic computations, etc.) 

8. Helpful hints concerned with the "best" way to do something (l.e., 
for linear elastic analysis use the option ‘constitutive to avoid 
unnecessary computations, etc.) 

9. What Information about the problem can be Inferred from other 
Information. All but the last type of knowledge appears In the 
user's manual. 

Of course, many of the first eight rule-types have been developed from 
talking with experts on the NESSUS code because the manual does not always 
provide all of the Information necessary to run the code. However, It does 
provide an excellent place to start. 

The overall design of the user Interface maintains In 
spirit, anyway, the three step process used by the NESSUS code for developing 
a data deck for the FEM processor Inputting the parameter data, the model 
data, and the Incremental data, if needed. Input to the pre-processor Is 
handled as a separate option In NESSUS/EXPERT. However, Inputting the 
parameter data Is not done Immediately at the start of a session because many 
of its values can be Inferred from the model data. Thus, the model data Is 
Input first, the necessary parameter values are determined by NESSUS/EXPERT 
and then the user Is given a chance to enter whatever other parameter data 
he/she deems necessary. Each of the three steps consists of the following 
rule-sets: 

o Rules to guide the questioning for required Information and to 
check Its correctness 

o Rules to handle the optional, keyword Input and to check Its 
correctness 

o Rules to handle a HELP facility 

o Rules to output the data to the screen so that the user can verify 
the data 

o Rules to check the completeness and consistency of the provided 
data 

o Rules to write the data to a file In the proper format 

Each of these groups of rules will constitute a separate 
portion of the knowledge base that we will refer to as rule-sets. They will 
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contain the various types of knowledge discussed In the previous section. 

Overall, NESSUS/EXPERT for NESSUS/FEM can currently be 
characterized as a menu-driven consultant. Some Input may come from 
previously prepared files, while other Input can be supplied by the user 
Interactively at the terminal. The basic tasks accomplished during each of 
the major phases In the system are described below. A block diagram of the 
system corresponding to this description appears In Fig. 2.25. 

1. Identify Problem : During this beginning step, the user Is asked to 
specify the name of the output file to be created and the type of 
data deck to be created. This Information Is then used by 
NESSUS/EXPERT to determine what should be done next. 

2. Define a Preprocessor Data Deck : If the type of data deck to be 
created is a pre-processor data deck then the system follows the left 
branch of the flow diagram In Fig. 2.25. Currently, NESSUS/EXPERT 

Is set-up simply so that such Information can be entered through a 
dialogue guided by the expert system so that everything that Is 
needed Is entered. Each data set must consist of five categories of 
Information: 1) RANDOM, 2) SELECT, 3) POINTS, 4) MEANS, and 5) 
DEVIATIONS. NESSUS/EXPERT simply prompts the user to enter all of 
this Information during the dialogue. The structure for consistency 
checking of the data before It Is written to the file Is available, 
but currently no rules have been Implemented. Any number of pre- 
processor data decks can be created during a given session. At the 
end of the session, the data Is written to the file specified 
Initially so that it can then be used by the pre-processor. 

3. Initial Dialogue : If, during the Initial Identification dialogue 
the user specified that a FEM data deck Is to be created, then the 
right branch of the system flow diagram given In Fig. 2.25 Is 
followed. The user is asked to provide some Introductory Information 
and to complete the minimum subset of model data categories which 
constitutes a valid data deck. This Information Is extracted through 
an initial dlalgue with the user which at the moment Is an unvarying 
sequence of questions for which the user must supply answers. This 
area of the code will eventually need significant expansion from the 
AI point-of-vlew. Currently It only contains a minimum amount of 
knowledge that was derived from the MHOST manual. Eventually, It 
will include more detailed expert knowledge that wll help to 
determine the categories that should be Included in this minimal data 
set based on some general questioning of the user. 

4. Input Model Oata : Most of the topics for the model data section of 
the NESSUS data deck are selected by having the user specify a topic 
by number or name from a large list of available topics. These 
topics correspond to the keywords used in the NESSUS code. Once a 
topic Is selected, NESSUS/EXPERT guides the user in Inputting the 
required Information associated with that topic either by hand of 
from an existing file. When Input Is completed, control Is returned 
to the main model data menu. Respecification and alteration of data 
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already provided Is possible at all times. Also available Is a help 
file on any of the topics that can be selected. Currently the 
knowledge that Is contained In this sefctlon has come mostly from the 
user manual for HHOST. However, when ambiguity or Inconsistency has 
appeared Information has been acquired directly from the coders and 
testers of HHOST. 

5. Input Random Variable Definitions and Perturbations : If a 
probabilistic FEM data deck Is being prepared In the HESSUS/EXPERT 
session, the user will be asked to provide a set of random variables 
and perturbations once he/she has completed the model data section of 
the data deck. In overall style, the data entry in this section Is 
handled In a manner highly consistent with previous sections of 
NESSUS/EXPERT. A certain set of keywords are needed, along with * 
their corresponding piece of data. The system provides guidance In 
filling in the values associated with the the required keywords 
either manually or from a file. First, the definitions are simply 
asked for, then input of the perturbation Information is guided by a 
parameter menu just like the model data section. As with all other 
sections of NESSUS/EXPERT, information can be corrected, deleted, or 
respecified at any time. This section currently embodies only the 
knowledge provided by Supplement II of the HHOST User's Hanual. 
However, this section will require much more attention In terms of 
providing support to the user In the form of an Intelligent aid for 
handling probabilistic geometric data In the coming year. 

6. Consistency Checking and Validation of the Data Deck : Consistency 
checking of the completed data deck is one of the more Important 
functions of NESSUS/EXPERT for it Is here that. much of the expert 
knowledge on how NESSUS works would be used to ensure a correct data 
deck. The goal of consistency checking Is to determine whether the 
information In the completed data deck Is consistent among all of the 
various categories. The rules encoded so far In NESSUS/EXPERT are, 
for the most part though sometimes very subtlely, contained in the 
MHOST User's Manual. Much of the knowledge has required clarification 
from either experts at SwRI or the original coders of the NESSUS 
system. When a problem Is detected In the Information provided in 
the data deck, the user is given a number of options for solving It, 
depending on the problem itself. Due to the power provided by a tool 
such as 0PS5, all errors will be detected In a very straightforward 
manner and If another Inconsistency Is created by fixing a problem, 
this is detected as well. The knowledge encoded in the system so far 
has emphasized compatibility between the parameter and model data, 
between the BFGS and ITERATIONS data, between the CONSTITUTIVE and 

the WORKHARD data, between the random variable data for a particular 
topic of the model data and that model data topic, between the 
perturbation and random variable data, and within the WORKHARD data 
itself. This section will continue to be expanded for the duration 
of the project as this is where much of the Intelligence of the 
NESSUS/EXPERT will reside. 

7. Output Data Deck : Once the data deck has been completed and verified 
as being consistent (to the extent that is currently possible by 
NESSUS/EXPERT), the data deck is printed out to a file. It is done 
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in the following order: 1) the header records and deck title card, 

2) the parameter data section, 3) the model data section, 4) the 
random variable section (If needed), and 5) the perturbations (If 
needed). The various sections are printed out In a suitable order 
(alphabetically or numerically as appropriate). This output goes to 
the file specified at the beginning of the session. Most of the 
basic structure of NESSUS/EXPERT exists now. What Is left to do In 
many cases Is to fill In the knowledge bases so that the coverage of 
NESSUS/EXPERT is complete. Other major additions left to be done are 
addressed In Section 3.1.3, Current Efforts on NESSUS/EXPERT. 

2. 1.3. 4 Rule Structure 

A production rule encodes knowledge about a problem In the 
form of IF-THEN statements also known as condition/action pairs. These 
production rules manipulate a set of data structures called objects. There 
can be an arbitrary number of these objects and each has associated with It a 
set of attributes and potential values for those attributes. 

The generic form of an 0PS5 production rule looks like 

the following: 

(p ex-rule (objectl attributel valuel attr1bute2 nil) (object2 
attribute3 <> value3) — > (make object4 attr1bute4 value4) 

(modify 1 attribute2 value2) ) 

The letter "p" just Inside the left parenthesis Indicates 
the beginning of the production rule. The rule's name Is "ex-rule". This 
allows the system to distinguish it uniquely from all other rules In the 
knoweldge base. The rest of the rule that occurs before the symbol " — >" 

Is called the left-hand-side (LHS) of the rule. It contains two 
conditions. The first is that there exist an objectl with an attributel of 
value valuel and an attr1bute2 with no value. The second is that there 
exist an object2 whose value for attribute3 Is not equal to value3. The 
portion of the rule following the symbol Is called the 

right-hand-side (RHS) of the rule. It contains two actions. The first 
creates a new object, called object4 with attr1bute4 of value value4. The 
second modifies the first object listed In the LHS of the rule (objectl) so 
that Its attr1bute2 has value value2. Thus, if this rule were to become 
true. It would result in modifying the world of objects and attributes in 
that specified way. 

In 0PS5 such rules are used during processing by a 
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technique called forward chaining. This means that the rules are data- 
drtwen. Data In the world (l.e., the objects and their attributes) change. 
Sudh changes cause some of the LHS's of the rules In the knowledge base to 
became true. One rule from this set of true rules Is selected through a 
metfcod called "conflict resolution" as the appropriate rule to activate, or 
fire. Firing causes the actions on the RHS of the rule to be performed 
resulting In changes to the data In the world making a different set of 
production rules In the knowledge base true. This process of forward chaining 
continues until Information Is needed from the user or no more production . 
rules can become true. If Information Is needed from the user, then this new 
Information can modify the data In the world, thus resulting In continuing 
the chaining process. If no more rules are true, then processing stops. 

One can represent fairly directly In 0PS5 the knowledge 
needed for NESSUS/EXPERT, such as information concerning a certain piece of 
parameter data for NESSUS. For example, the object could simply be called 
parameter-data. Its attributes could Include Its name. Its parameter-value 
names, and related model data names. The parameter data name's value could 
be ELEMENTS, its first parameter-value (element type) could be 7, and the 
related model data names would Include ELEMENTS. A rule could then be 
devised that, based on the fact that the parameter data's name Is ELEMENTS 
and Its first parameter value is 7, can determine which pieces of model 
data are needed to run the problem correctly. The rule might look 
something like the following when converted Into English: "IF there Is an 
object called parameter-data, whose name is ELEMENTS and whose first 
parameter-value is 7, THEN the model data whose name is ELEMENT Is also 
needed. This Is a fairly obvious and simple rule, but they can become very 
complex, depending on what knowledge must be represented. The result of 
thfs rule Is that If parameter data called ELEMENT exists In the data deck, 
tbetn the corresponding model data called ELEMENT must also exist. This Is 
a simple example of how consistency checking of the data deck can be done 
ustnig 0PS5 rules. 

2.1.4 Verification Studies 

2.1.4. 1 Objectives of Verification Efforts 

The basic objective of the verification effort Is to 
apply the methods developed and implemented In NESSUS family of computer 
programs to the analysis of actual space propulsion system components. The 
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typical components to which the methods will be applied Include a turbine 
blade, a high pressure duct, a lox post and a transfer tube duct liner. The 
verification efforts would cover a wide range of analysis options developed 
and Implemented In NESSUS codes. 

The knowledge gained In the verification efforts will be 
Implemented In NESSUS expert system. The verification effort Is broadly 
divided Into simple verification and complex component verification 
analysis. Since NESSUS Is In a state of continuous development during the 
contract, the simple verification studies are designed to meet the following 
objectives. 

The simple verification models exercise the element 
types, the typical random variables, the range of perturbation of each random 
variable and various solution strategies for a particular component but on a 
simplistic model. These studies differ from validation studies by the fact 
that they are specifically targeted for each component analysis. 

The results of the simple verification studies aid in 
establishing confidence In the code, Identify Its limitations In user 
Interface, as well as analysis capabilities when applied to analysis of 
practical components. They also result in correcting element deficiencies and 
devise solution strategies that will be effective when analyzing full scale 
verification problems. The full scale verification problems on the other 
hand. If possible, are conducted on existing production finite element models 
and are typically expected to be much more computationally intensive requiring 
large main frame computing facility. 

2. 1.4. 2 Scope of Verification Efforts 

The space propulsion system components are subjected to 
environments with many random variables. Due to the difficulties in the 
Instrumentation of high energy, high pressure and temperature systems, many 
variables are not wel 1-characterlzed. Nevertheless, many components are 
subjected to severe environments. The current design philosophy Is to analyze 
and design the components based on worst conditions using state-of-the-art 
deterministic analysis methods. The environments and conditions under which 
many space propulsion system components operate lead to structural analysis in 
the non-linear analysis domain. These structural analysis non-linearities can 
be due to material property or due to geometric changes or due to contact 
boundary conditions. Detailed discussion of the environments and 
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deterministic analysis techniques to which some typical space propulsion 
system components such as turbine blades, lox posts, transfer tubes, high 
pressure ducts, nozzle feed lines, and main combustion chamber walls are 
subjected to were described In detail In the first year report. 

The composite loads spectra contract and probabilistic 
structural analysis contract are bold and challenging attempts to extend 
advanced deterministic structural analysis methodologies Into probabilistic 
structural analysis domain. Developments under the PSAM contract are 
Implemented Incrementally Into the HESSUS program during the five year 
contract period with Increasing levels of analysis sophistication each year. 
Due to scheduling constraints, all analysis options available In NESSUS can 
not be applied to every component. Thus, a strategy has been developed In 
which the component, the type of structural analysis, random variables and the 
area of emphasis are chosen to be consistent with code development. This has 
been achieved In a probabilistic structural analysis domain for each component 
consistent with primary deterministic analysis requirement for each 
component. The scope of the verification studies achieves these objectives 
for each component In the order listed below: 

1. Turbine blade 

2. High pressure discharge duct 

3. Lox post 

4. Transfer duct 

2. 1.4.3 Turbine Blade Component Random Variables 

The high pressure fuel pump turbine blade has been c! osen 
as the first component to be analyzed by NESSUS finite element code. The 
analysis options and random variables chosen are consistent with the state of 
program development. The random variables that will be exercised on turbine 
blade analysis are: 

1. Material property variations and orientations. 

2. Geometry changes. 

3. Centrifugal loads. 

Pressure loads. 


4. 



79 


5. Temperature loads. 

The strategy of the treatment of the random variables are first presented 
followed by the results of the simple verification studies. 

2. 1.4.4 Material Property Variations 

The most commonly used turbine blade materials are 
nickel-based super alloys. Directionally solidified Mar-M-246 (Hf) Is used In 
space shuttle main engine high pressure turbopump turbines. The single 
crystal PW1480 material Is being evaluated for future use In the engine. 

These materials are anisotropic In nature and exhibit strong directionally 
oriented properties. As an example, for the PW1480 material at room 
temperature, the elastic modulus In the 111 plane can be as much as 250% 
greater than the modulus In 001 plane (Fig. 2.26). Thus, any perturbation of 
material orientation affects the blade stiffness and thereby Its static and 
dynamic response. The material orientation angle Is one of the random 
variables chosen In probabilistic structural analysis of turbine blades. 
Treatment of material orientation angle In single crystal blades Is easier 
when compared to Directionally Solidified (DS) material blades. This Is 
because the DS blade material typically contain a random number of crystals In 
each blade, (usually from 3 to 10), the volume of which can be random, with 
each crystal having its own material axis orientations. The single crystal 
materials, on the other hand, contain only one crystal but the orientation 
angles can vary slightly along the length of the blade based on crystal growth 
direction. 

A typical statistical data of the distribution of the 
primary material axis orientation to the stacking axis from a set of hundred 
blades as measured at the base of firtree Is shown in Fig. 2.27. Statistical 
analysis of data Indicates a normal cumulative distribution provides a 
reasonable good fit of data. However, since blades having a cone angle of 
greater than 10° are rejected, the cumulative distribution function for the 
accepted blades is a truncated one modified as shown In Fig. 2.27. 

Perturbation of material orientation angles Is achieved 
In NESSUS by designating the orientation angles as a random variable. The 
studies of the perturbation of material orientation angles and the behavior of 
the numerical algorithm Is discussed later In the section. 

The other factor that might be considered in the material 
property variations is the scatter in elastic constants themselves from 
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Fig. 2.26 Typical Elastic Modulus Variation From 001 Plane to 111 Plane for 
Nickel-Based Super Alloys 
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specimen to specimen at a given temperature. In general. Insufficient amount 
of elastic properties data points exist at each temperature to do a good 
statistical analysis to accurately characterize the variations. However, this 
variation for elastic constants at a given temperature Is small. A typical 
set of elastic constants for wide range of temperatures for PW1480 material Is 
shown In Table 2.6. 


Table 2.6 

Elastic Constants for PW1480 Material as a Function of Temperature 



-400°F 

70°F 

1400° F 

2000 °F 

E 

19.96E6 

18.38E6 

14.75E6 

11.0E6 

G 

20.50E6 

18.6366 

15.2786 

12.82E6 

n 

0.376 

0.386 

0.395 

0.416 


The material property for anisotropic material Is currently Input to the code 
explicitly by specifying completing the material 1 D matrix (s=De). However, 
for PW1480 material in the principal material orientations, a set elastic 
constants that can completely characterize the elastic response can be 
specified by E, n, and G. Thus, new features will be added to the code for 
specifying these constants (Instead of the full D-matrlx) and perturbations of 
them to calculate the response due to material property variations. The 
option of perturbing each coefficient of the full D-matrlx Is postponed to 
later releases of the code. The Issue of building In rules In the NESSUS 
expert system to avoid material property perturbations that violate the laws 
of physics such as non-positive definiteness of the matrix will be addressed. 

2. 1.4.5 Geometry Changes 

Because of the criticality of the component, every 
turbine blade that Is used In an engine is subjected to quality inspection 
procedure for adherence to the design geometry. The blades that are used In 
space propulsion systems are typically short and compact, 0.5" to 3.0" In 
length when compared to turbine blades used In air breathing engines. The 
specified tolerance is a function of the manufacturing method. For cast 
blades, the tolerances are usually of the order of 0.005". Many turbine 
blades, including the kind used in the Space Shuttle Main Engine, are of cast 
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type with machined flrtree which forms the mechanical attachment to the 
disk. The measured geometrical variations found In these blades generally 
fall In the category of relative twist of the blade (Fig. 2.28) and lateral 
shift of the profile within the tolerance envelope. This Is presented In the 
form of center of gravity shift (x, y, coordinates) for a set of about seventy 
blades (Fig. 2.29). An analysis of the measured data Indicates that a 
majority of geometrical variations from blade to blade occur when the flrtrees 
are machined. The net effect of geometric variations Introduced In this 
machining step Is a rigid body shift of the airfoil, shank and platform 
relative to the stacking axis which runs at the center of flrtree. Thus, the 
strategy that will be adopted for the perturbation of geometrical quantities 
for turbine blades will be the perturbation of nodel coordinates of the finite 
element model resulting from rigid body rotation about x, y and z axes 
rotations. 

2. 1.4.6 Centrifugal Load Variations 

Centrifugal load Is one of the primary loads on turbine 
blades. It contributes to a major share towards the mean stress level and 
thus plays a critical role in fatigue life calculations. The centrifugal load 
varies as the square of the turbine speed. The speed profile of high pressure 
fuel turbopump In SSME is shown In Fig. 2.30 which closely follows the engine 
thrust profile. An expanded trace of measured speed between 32000 to 36000 
rpm from a pump signature test Is shown In Fig. 2.31. Here, the power level 
was reduced 1* per three seconds of test. 

Random speed oscillations can be seen about a mean from 
this data. Detailed study of test to test variations furnishes a good 
statistical database for this data. It Is a level 1 type of probabilistic 
loading In that randomness of centrifugal load Is spaclally homogeneous for 
the finite element model. The engine balance models Indicate that 2s speed 
variations at steady state power level for the SSME fuel pump Is about 400 rpm 
out of 36600 rpm assuming a normal distribution. It Is planned to use the 
actual processed test data from engine tests for the probabilistic structural 
analysis. The benefit of the results from the composite load spectra 
development contract will be utilized for all loads subject to their 
availability. 
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Fig. 2.31 Expanded Trace of Speed Profile From a Pump Signature Test 
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2. 1.4.7 Steady- State Pressure Loads 

The steady state pressure loads on turbine blades Is a 
function of flow conditions at Inlet and outlet of the turbine. Detailed 
measurements of turbine blade surface pressures and temperatures from actual 
engine tests are unavailable. There are a number of measurements such as 
preburner chamber pressure, downstream turbine discharge pressure and 
temperature (downstream of turnaround duct) and pump head raise 
measurements. There are a few measurements from instrumented turbopumps for 
temperature In the stators (nozzles) upstream of turbine blades. Thus, the 
fluctuation of static differential pressure on the turbine blade between 
pressure and suction faces will be a calculated quantity obtained from 
Indirect measurements and theoretical engine models. 

The type of stochastic modeling of pressure load on a 
turbine blade Is closely related to the design features of the turbine. For 
the chosen high pressure fuel turbopump component, the design features are 
Illustrated In Fig. 2.32. A notable feature Is that this turbine has a 
secondary flow circuit for cooling the rotating hardware and includes cooling 
of the shank portion of the turbine blades. ..This cooling circuit affects the 

i 

pressure in shank portion of the turbine blade. Thus, the pressure load on 
turbine blade will be treated as a random field. Level II type modeling. It 
Is planned that the statistics of the differential pressure variation for the 
airfoil will be correlated through turbine torque variation. The shank 
pressure variations will be correlated to coolant pressure variations. 

Typically, the pressure Information will be availabTe at 
three or four streamlines or cross sections which will be Independent of the 
particular finite element model. The pressure at model node locations for a 
particular model will then have to be obtained through interpolator codes. 

2. 1.4.8 Blade Temperature Loads 

The temperature loads plays a critical role In turbine 
blade analysis. For space propulsion systems of L0X/LH2 systems with staged 
combustion process, the range of temperatures can be very high in a duty 
cycle. For example. In SSME during one mission duty cycle, the blades will be 
a temperature range from 2200/R to 200/R. While it is virtually impossible 
to measure turbine blade temperatures in an actual engine, first stage stator 
(nozzle) temperature data from a few instrumented turbopumps is available. 
While temperature transients cause the worst case stresses when compared to 
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steady-state, the Initial scope of the probabilistic structural analysis of 
turbine blade will be limited to steady-state loads. For this reason, the 
discussion Is limited to solution strategies for the treatment of 
probabilistic temperature loads at steady-state. 

Just as In the pressure case, the characteristics of the 
temperature random variable is closely related to the design features of the 
turbopump. For a high pressure fuel turbopump, the coolant flows around the 
shank. In actuality, the coolant and hot gas flows around the shank are very 
complicated. The hardware shows large variations In oxidation discoloration 
(which Is a rough Indication of temperature) from pump to pump. Indicating 
that as the various seals wear they affect the flow circuit resistances and 
thereby temperature in the shank region. Thus, the developed probabilistic 
structural analysis methodology should be able to handle large local 
perturbations In temperature. On the other hand, the airfoil temperatures at 
steady state is essentially the hot gas temperature. Typically the shank area 
has a large thermal gradient when compared to the airfoil as shown In Fig. 

2.33 and Fig. 2.34. The platform of the turbine blade itself Is nearly 
Isothermal at steady-state. Thus for the probabilistic structural analysis of 
HPFTP turbine blade, the temperature will be treated as a random field with 
varying statistical characteristics in airfoil, platform and shank. Thus, 
stochastic modeling of temperature Is a Level II type modeling. 

2. 1.4.9 Deterministic Verification Solutions 

Simple models, Fig. 2.35 through Fig. 2.37 comprised only 
of solid elements were exercised in NESSUS/FEM to understand and verify the 
performance of basic solid element as Implemented In NESSUS. Several random 
variables were also exercised with typical range of perturbations that will be 
used In component verification studies. First, the deterministic results 
obtained from NESSUS are discussed, followed by perturbation analysis 
results. All the exercises were conducted on an anisotropic beam 
representative of PW1480 material properties at room temperature. 

Considering centrifugal load first, model shown In Fig. 
2.35 was exercised for both with one of the model axis as the axis of rotation 
as well as off axis rotation for hinged condition. The program results 
exactly match the theoretical calculated radial loads. 
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Fig. 2.33 Typical Temperature Gradient In the Shank Region for High Pressure 
Fuel Turbine Blade 
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Fig. 2.34 Typical Temperature Gradient In the Airfoil Region for High 
Pressure Fuel Turbopump 
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Fig. 2.35 Solid Element Beam Model A, 2 x 20 Elements 



Fig. 2.36 Solid Element Beam Model B, 4 x 20 Elements 
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Fig. 2.37 Solid Element Beam Model C, 8 x 40 Elements 
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Considering the pressure loads next, the models In Fig. 
2.35 - Fig. 2.37 were exercised for uniformly distributed pressure load and 
constant moment condition. The results of the uniformly distributed pressure 
load Is presented In Table 2.7. 


Table 2.7 

Uniform Pressure Loading Results Cantilever Beam 
FEM Results/Theory Ratio 



Simple Beam 
Theory 

Model A 

Model B 

Model C 

Tip Deflection 

1.0 

0.69 

0.74 

0.877 

Fixed End Stresses 

1.0 

0.51 

0.66 

0.867 


The basic solid elements as Implemented currently In 
NESSUS Is a strict elght-noded Isoparametric element. It Is known that these 
elements are stiff when they encounter pure bending situations and require a 
fine mesh to obtain good results. There are several approaches possible to 
Improve the performance of this element. One of the well-known approaches is 
the Introduction of additional modes such as (1-r 2 ), (1-s 2 ), (1-t 2 ) for the 8- 
noded brick elements. While the introduction of these functions Improves the 
performance dramatically for pure bending cases, they also violate 
compatibility and do not pass the patch test for arbitrary shaped 
quadrilaterals. Further, the performance deteriorates for arbitrary 
quadrilaterals. The problem of the patch test failure was subsequently cured 
by evaluating the contribution of the incompatible modes to the jocobian 
matrix at the centroid. It has been found that the resulting element gives 
superior performance to the original Incompatible element. 

The other approach to make the element flexible Is 
through the use of reduced Integration quadrature. The two concepts that are 
used are fully or uniformly reduced quadrature and selective reduced 
Integration quadrature. Recent studies demonstrating the equivalence of a 
class of mixed models with reduced/selective Integrated elements In linear 
elasticity has elevated the reduced integration approach from "tricks" to 
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legitimate methodology. However, the Important considerations In the use of 
the methods are the Insufficient rank of the matrix In the fully reduced 
method and the extension of the methodology to anisotropic cases In selective 
reduced method. The fully reduced quadrature Is available In HESSUS without 
the hour-glass control. The fully reduced quadrature results In spurious 
modes, and therefore must be used with caution. For static analysis, 
computations using fully reduced Integration scheme may be possible depending 
upon the boundary conditions providing stability to the problem. However, for 
transient dynamic analysis, hour-glass viscosity control to suppress the 
spurious modes Is a necessity to obtain accurate results. 

One of the principal deficiencies of the selective 
Integration procedure or recently the B approach as normally Implemented Is 
that it is limited to Isotropic case. For turbine blade applications, the 
material Is anisotropic and the D-matrlx Is fully populated for general 
material orientation. The use of standard selective reduced Integration 
schemes to anisotropic cases Is ambiguous. Thus, It Is desirable to Implement 
extensions to selective Integration schemes or to the B approach in the 
context of displacement formulation to cover anisotropic cases. The 
additional benefit of such a procedure would be Its extension to nonlinear 
problems where tangent modulll always exhibit anisotropic character. Several 
temperature gradient solutions were also conducted on models Fig. 2.35 through 
Fig. 2.37 for the anisotropic material element. One of the notable features 
of the PW1480 material is that while its elastic properties exhibit strong 
directionally dependent properties, the coefficient of thermal expansion is 
nearly Isotropic. The results of the temperature solution are presented In 
Table 2.8. 
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Table 2.8 

Temperature Gradient Solution 1000° Through Thickness 



Simple Beam 
Theory 

Model A 

Model B 

Model C 

Tip Deflection 

(ratios) 

1.0 

1.17 

1.02 

1.12 

Stress 

(absolute 

values) 

0 

128000 

68000 

33786 


The maximum perturbations from 001 to 111 and vice-versa, 
were tested to check the convergence characteristics under maximum elastic 
property changes resulting from material orientation.. In practice, material 
orientations are not allowed to differ more than ±10° from the primary 
direction. Thus, the Newton-Raphson method Is expected to be adequate for 
material orientation perturbations for component verification. The same 
strategy should also be adequate for material property variations also as they 
are typically very small for single crystal blades. 

Perturbation studies on geometrical changes are next 
addressed. The rigid body rotation about the base of the cantilever type 
geometric variations found In SSME turbine blades were earlier discussed. The 
greatest effect of this type of variation Is In the contribution due the 
centrifugal load to the stresses due to change In the center of mass 
location. Two studies were conducted on the Model A ( Fig. 2.35) where the 
geometrical perturbations were 1 degree and 10 degree rotational shift about 
the base of the rotating beam, for 1 degree perturbation, the default Newton- 
Raphson method converges for normal englneerlngllmlt of acceptable residual 
load errors, however, when the residual load vector Is tightened to the order 
of IE-5 of the total centrifugal load, the Newton-Raphson technique exhibits 
convergence and then divergence characteristics. However, when sealant 
Iteration option Is used, the algorithm exhibits uniform convergence and 
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converges to the tight tolerances (IE-5 of the total load) in three 
Iterations. In actual turbine blades, variations of less than 1 degree tilt 
are expected. Thus, the available solution strategy appears adequate to 
handle geometric perturbations. 

Due to the convergence behavior of Newton-Raphson 
technique for 1 degree perturbation, another case with a 10 degree 
perturbation was run. This case, the Newton-Raphson technique diverges from 
the start. The results are still under study. One fo the features of the 
test problem Is the state of stress and centrifugal load In body fixed 
reference frame do not change due to perturbation. However, the global 
location of the body Is different when measured from determinate reference 
frame after perturbation. The question of how large a perturbation the 
Implemented solution strategies can tolerate will be studied further. 

At the current state of development, NESSUS/FEH Is 
applicable for linear analysis only. Thus, the perturbation of loads such as 
centrifugal and ressure loads amount to resolving the linear problem for a new 
load case with the old stiffness matrix. Irrespective of the magnitude of the 
perturbation of centrifugal and pressures, solutions converged In test cases 
In two Iterations using Newton-Raphson method. Perturbation of loads end 
convergence have a greater bearing in the nonlinear analysis. The simple 
verification studies will continue to Improve element and algorithm 
performances under a variety of conditions. Some of the Improvements under 
development from the verification studies are described In the current efforts 
chapter of NESSUS/FEM. The results of the study will be used In component 
verification analysis of the turbine blade. 

2.1.4.10 Perturbation Verification Studies 

Perturbation verification studies were conducted on the 
model shown In Fig. 2.35. The random variables exercised to date Include: 

i 

1. Material orientation angle 

2. Nodel coordinates 

3. Pressure 


4. 


Centrifugal Load 
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The perturbation algorithm relies on established 
predictor-corrector methods used In nonlinear finite element analysis. There 
Is no one Iterative method that exists that exhibit Ideal convergence 
characteristics as well as be cost effective In all situations. The solution 
strategy to be used Is a function of the type of nonlinearity at hand. The 
methods that have been developed for nonlinear finite element analysis Include 
full Newton, Quasi Newton, and Newton Raphson techniques. All the above 
techniques are available In NESSUS at a global level common to all 
perturbations within a run. 

The logic for choosing the solution strategy should 
primarily depend on the rate of convergence and cost of the solution. A 

necessary condition for convergence for all the Iterative methods Is the exact 

calculation of residual load vector at each Iteration. They all differ in the 
evaluation of predictor, the trial stiffness matrix used. In full Newton, the 
tangent stiffness is evaluated at every Iteration. In the modified Newton- 
Raphson, the original stiffness matrix or the matrix at the start of the 
Increment Is used. In Quasi-Newton methods, the stiffness matrix Is updated, 
but numerical strategies are used to reduce the amount of computations (update 
of stiffness matrix without inversion) than It would be If a full Newton 

method (requiring a full matrix Inversion) was used. The Initial exercises In 

the perturbation examples use the default Newton-Raphson method In the code. 
Other solution strategies were used only when divergence was encountered while 
using Newton-Raphson method. 

The material angle perturbations are first addressed. 

The model (Fig. 2.35) was exercised for material axis variations In the 
presence of pure axial load. The objective of the studies were to test the 
convergence characteristics. One of the considerations was the study of the 
performance of the default Newton-Raphson method under perturbations that 
stiffen the structure. This can happen In turbine blade analysis when 
material orientation variations can result In a stlffer blade In the primary 
radial direction. 

The study exercised the model In Fig. 2.34 with 
perturbations about the deterministic state resulting In stlffer or softer 
structure with varying magnitude. The results are summarized in Table 2.9. 
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Table 2.9 

Material Orientation Angle Perturbation 
Axial Load Results 


No. of Iterations 
for Residual Load 
Environment 

Deterministic Amount of Perturbation % of Applied Load 

State About Deterministic State Convergence 1 0.1 0.01 


001 

+ 10° 

yes 

4 

8 

16 

001 

To match 111 plane 
(36° + 45°) 

no 

- 

- 

- 

111 

+ 

o 

o 

yes 

2 

3 

7 

111 

To match 001 plane 
(36° + 45°) 

no 

- 

- 

- 


The maximum perturbations from 001 to 111 and vice-versa, 
were tested to check the convergence characteristics under maximum elastic 
property changes resulting from material orientation. In practice, material 
orientations are not allowed to differ more than +10° from the primary 
direction. Thus, the Newton-Raphson method is expected to be adequate for 
material orientation perturbations for component verification. The same 
strategy should also be adequate for material property variations also as they 
are typically very small for single crystal blades. 

Perturbation studies on geometrical changes are next 
addressed. The rigid body rotation about the base of the cantilever type 
geometric variations found in SSME turbine blades were earlier discussed. The 
greatest effect of this type of variation Is In the contribution due the 
centrifugal load to the stresses due to change In the center of mass 
location. Two studies were conducted on the Model A (Fig. 2.34) where the 
geometrical perturbations were 1 degree and 10 degree rotational shift about 
the base of the rotating beam. For 1 degree perturbation, the default Newton- 
Raphson method converges for normal engineering limit of acceptable residual 
load errors. However, when the residual load vector is tightened to the order 
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of IE-5 of the total centrifugal load, the Newton-Raphson technique exhibits 
convergence and then divergence characteristics. However, when sealant 
Iteration option Is used, the algorithm exhibits uniform convergence and 
converges to the tight tolerances (IE-5 of the total load) In three 
Iterations. In actual turbine blades, variations of less than 1 degree tilt 
are expected. Thus, the available solution strategy appears adequate to 
handle geometric perturbations. 

Due to the convergence behavior of Newton-Raphson 
technique for 1 degree perturbation, another case with a 10 degree 
perturbation was run. This case, the Newton-Raphson technique diverges from 
the start. The results are still under study. One of the features of the 
test problem Is the state of stress and centrifugal load In body fixed 
reference frame do not change due to perturbation. However, the global 
location of the body Is different when measured from determinate reference 
frame after perturbation. The question of how large a perturbation the 
Implemented solution strategies can tolerate will be studies further. 

At the current state of development, NESSUS/FEM Is 
applicable for linear analysis only. Thus, the perturbation of loads such as 
centrifugal and pressure loads amount to resolving the linear problem for a 
new load case with the old stiffness matrix. Irrespective of the magnitude of 
the perturbation of centrifugal and pressures, solutions converged in test 
cases in two iterations using Newton-Raphson method. Perturbation of loads 
and convergence have a greater bearing in the nonlinear analysis. The simple 
verification studies will continue to improve element and algorithm 
performances under a variety of conditions. Some of the Improvements under 
development from the verification studies are described In the current 
elements chapter of NESSUS/FEM. The results of the study will be used In 
component verification analysis of the turbine blade. 



103 


3.0 CURRENT EFFORT 


3.1 NESSUS/FEM 

Two different approaches have been proposed for the extension of the 
NESSUS perturbation algorithms to Inelastic problems. The first approach 
calls for continuing the development within the displacement formulation 
used In the first year PFEM effort. Extension of the displacement formulation 
to Inelastic analysis in NESSUS/FEM will require a major reorganization of the . 
internal data structures within the code. The second approach calls for the 
adoption of a mixed Iterative formulation, which would preserve the Internal 
data structure of the present code. The development and Implementation of 
appropriate perturbation algorithms for Inelastic analysis will be started as 
soon as a decision is reached regarding the finite element formulation adopted 
for future PFEM development. 

The development of a finite deformation kinematics algorithm for 
NESSUS Is currently well underway. The adopted formulation utilizes an 
updated Lagranglan mesh description, with a constitutive relation based on 
the Green-Naghdl rate of Cauchy stress and rete^of deformation. Although 
the development of nonlinear displacement and strain modeling capability Is 
not required In NESSUS/FEM until FY88, MARC has taken advantage of the 
development of a similar capability for the MHOST code. The finite 
deformation algorithms being developed for the MHOST code will be added to 
the main development version of NESSUS/FEM In a very near future. 

An enhanced continuum-based plate/shell element with surface node 
definition Is currently under development at MARC. This element is 
envisioned as an eight-node brick with assumed strain modes based on the 
exact bending solution for an elastic Isotropic material. The approach Is 
expected to result In a non-locking element with enhanced bending behavior 
which can be distorted to a high aspect ratio (h/L - 1/10) In order to 
model moderately thick plate and shell-like structures. An early version 
of this element for use In linear elastostatics should be available In time 
for the 2/1/87 code delivery. 

A revised format to allow specification of surface pressures and edge 
tractions on a nodal basis will be developed and tested. Changes will be 
Implemented to allow the degeneration of continuum- type elements to form 
triangles, wedges and tetrahedra. This will require changes to the strain 
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smoothing procedure available in NESSUS 1.1. The new smoothing procedure 
will then be tested for robustness. 1 

The test cases proposed In the preliminary plan for validation of the 
NESSUS code are being exercised at MARC. In addition, MARC Is In the 
process of compiling a standard list of test problems that will be used to 
exercise all versions of NESSUS shipped from MARC. These problems range Tn 
size and complexity from small one element tests to Irregular element 
meshes of a few hundred degrees-of-freedom. 

3.2 NESSUS/FPI 

Testing of the new CDF estimation procedure (Section 2. 1.2. 7) and 
the validation of the NESSUS code Is In progress. The exact solutions of the 
validation test problems have been obtained for the first five problems 
(Section 2. 1.2.8). Solutions for the remaining problems will be obtained In 
the current year. These solutions will be used to compare results generated 
from the NESSUS/FPI. By using the perturbation data base generated by MARC 
(perturbation solutions are now available for the validation problems 1, 3 and 
5), the new CDF estimation procedure will be used to continue the validation 
of the NESSUS modules and the solution procedure. The solutions will require 
additional runs of the deterministic FEM solutions and. If necessary, 
additional perturbations. 

Effort In Integrating the NESSUS/FPI with the NESSUS/PRE, NESSUS/FEM and 
NESSUS/EXPERT Is In progress. The basic structure of the expert system code 
NESSUS/EXPERT Is in place and operational. The emphasis during the next year 
will be to make the code easier to use by the engineer. 

One of the difficulties Identified In conducting probabilistic structural 
analysis on systems with a large number of random variables is developing a 
method of efficiently entering the random variables Into the computer. For 
the analysts to enter a separate probabilistic data base would be time 
consuming and error prone. The approach being pursued Is to use the existing 
data base for the structural model along with the NESSUS/EXPERT to query the 
user as to which variables are random. Distributional Information and the 
degree of correlation will also be provided at this time. With this 
information, NESSUS/EXPERT can generate an input file for the FORTRAN code 
NESSUS/PRE. 

The user will now have to exit NESSUS/EXPERT to run NESSUS/PRE. 

However, prior to exiting, NESSUS/EXPERT will save a data file of the 



105 


model and the random variables. NESSUS/PRE transforms a set of 
correlated random variables to a set of uncorrelated random variables 
using the eigenvector transformation method. This set of uncorrelated 
variables are saved In a file. Finally, the user will have to enter 
HESSUS/EXPERT again. NESSUS/EXPERT will retrieve the previously stored 
files and generate a complete NESSUS/FEH file which Includes the 
structural model data, the random variables data and the perturbation 
settings. 

3.3 NESSUS/EXPERT 

Now that the basic structure and approach to NESUSS/EXPERT has been 
designed and Implemented, emphasis Is turning to an evaluation of this' Initial 
prototype to determine what Is good and bad about It. Work will also proceed 
on extending the knowledge base to Include knowledge of all keywords listed In 
the MHOST User's Manual. Finally, once the results of the prototype * ~ 

evaluation are completed and Implemented, work will begin on handling the 
probabilistic data In a more natural and Intelligent manner. 

Extensive discussions between the experts on the use of NESSUS and the 
knowledge engineer Implementing NESSUS/EXPERT have already begun, results so 
far Indicate that some changes to the basic control structure need to be made 
In order to take advantage of some overlap In the use of certain data In 
different sections of the input data deck. The result will probably be a 
major change to the overall flow diagram given In Figure 2.25. However, 
because of the use of a very high level language such as 0PS5, the necessary 
changes should not be difficult to make. 

Work on enhancing the knowledge base will not proceed until the basic 
changes to the prototype flow of control have been made. The major source of 
knowledge will be the MHOST User's Manual and the knowledge will emphasize the 
sue of keywords. More held files and consistency checking rules will also be 
added as the project progresses. When the Information In the User's Manual Is 
Incomplete or ambiguous, knowledge will be solicited from human experts on the 
use of MHOST. 

Handling the input of the probabilistic data In a natural and Intelligent 
way will require some research on what the best Interface might be. 

Currently, the method of Inputting of the data Is the same as for the model 
data. The knowledge that this section of NESSUS/EXPERT contains is simply 
Information about the keywords pulled mainly from the User's Manual. However, 


this method Is not very helpful or efficient for the user to work with when 
entering such Information. Possible enhancements Include providing some 
graphic aids that can Illustrate various permutations on an element and some 
Intelligence of probability as It relates to FEM so that HESSUS/EXPERT can 
make many of the decisions and perform many of the needed calculations Itself, 
rather than making the user do them. 
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Introduction 

An algorithm for the efficient computation of the elastostatic 
response of a perturbed system discretized vith finite elements has been 
proposed [2] and implemented in the NESSUS code as part of the PSAM 
development effort. Although this algorithm has been successfully used in 
sensitivity studies of several structural systems vith random parameters, 
recent experience has indicated loss of stability for seemingly "small" 
perturbations in some problem classes. These problems typically involve 
approximate constraint equations which are embedded in the stiffness of 
tne unperturbed problem, and perturbations which result in the 

modification of these constraint equations. 

Finite element formulations for constrained problems using Lagrange 
multipliers and the penalty method have enjoyed widespread use in the 
recent past and have played an essential role in the development of 
successful methods for certain classes of problems. The literature on 
this subject is extensive and includes applications to: 

o The analysis of plate and shell structures allowing shear 

deformation [7, 12, 16]. 

o Incompressible elasticity, e.g., rubber-like materials [5, 8]. 

o Deviatoric rate-independent plasticity [10]. 

o Incompressible flows, e.g., Stokes and Navier-Stokes equations 
[5, 14, 15], etc. 

A fundamental assumption in the development of the perturbation 
algorithm in [2] is that the response of the unperturbed system 
constitutes a "good" approximation to the response of the perturbed 
system. This will not, in general, be the case, if the prescribed 

perturbation results in a noticeable change in the constraint equations 
present in the unperturbed system. Violation of this condition will often 
result in loss of stability and failure to converge. Thus, the presence 
of constraint equations in the finite element formulation may impose 
limits in the size of some perturbation parameters which are not 
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immediately obvious. Additional analysis must then be performed in order 
to determine exactly vhat constitutes a "small" perturbation. 

Perturbation methods based on Taylor series expansions about the 
unperturbed solution have also been proposed by several researchers [3, 4, 
11] and it is natural to ask hov these pathologies manifest themselves in 
the solutions obtained by these methods. It can be shown that the 
displacement correction obtained in the first iteration is identical to 
the first-order term in the Taylor series expansion, and the one obtained 
in the second iteration is identical to the second-order term in the 
series. Thus, the rate of convergence of the iterative algorithm is 
closely related to the errors resulting from truncation of the Taylor 
series. One advantage of the iterative algorithm in [2] is that an error 
estimate (the force residual) must be computed and is available at every 
step of the iteration. 

The Transverse Shear Constraint 


The classical Poisson-Kirchhof f theory i, of plates requires 
continuity of displacement, as does the classical Bernoulli-Euler beam 

theory. However, the development of compatible C* interpolations in 
multi-dimensional cases is not straightforward, and considerable efforts 
and ingenuity were invested in the development of the first generation of 
finite element formulations for thin plate and shell problems 1 13]. 

In recent years, the Reissner-Mindlin theory of plates, which can 
accommodate transverse shear strains, has enjoyed widespread use. In this 

formulation, only continuity of displacements is required, allowing the 
construction of far simpler and less restrictive interpolation schemes. 
As a result, finite element formulations for medium-thick plate and shell 
problems have been developed, which retain accuracy even for thin plate 
and shell situations (7, 12, 16]. However, as the thin limit is 
approached, the "pure bending" modes dominate the solution, resulting in 
the emergence of penalty constraint terms in the stiffness equations. 
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The fundamental aspects of the problem may be observed in the one- 
dimensional analog of the Reissner-Mindlin plate theory, i.e., the 
Timoshenko beam. Thus, the Timoshenko beam theory may be used as a 
simpler, more manageable model exhibiting the pathologies that afflict 
Reissner-Mindlin plate theory. The total potential energy, including 
shear deformation, for an elastic beam of rectangular cross-section with 
thickness t and width b may be written as 


1} 
2 J 
o 


Ebt' 

12 


'dx' 


dx 


i} 
2 J 
o 


tcGbt 


<£ - 0)2 


dx - J wqdx 
o 


( 1 ) 


In this form, the first integral corresponds to the "pure bending" energy 
in the beam, whereas the second integral represents the shear deformation 
energy, and the third and last term accounts for the work done by the 
applied transverse loading. As the thickness t is reduced, the bending 

stiffness (Ebt 3 /12) will decrease much faster than the shear stiffness 
(xGbt/2) . In the limit, the shear stiffness term becomes a Lagrange 
multiplier enforcing the condition that 


dv 

dx 


e 


which is the assumption made a priori in Bernoulli-Euler beam theory that 
the rotation is the derivative of the transverse displacement. 


The Discretized Problem 


The finite element formulation for the Timoshenko beam problem using 
linear interpolations for both translational and rotational degrees-of- 
freedom produces an element stiffness matrix of the form 
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mhwre h is the element length, E and G are the elastic and shear moduli, I 
and A are the cross-sectional moment of inertia and area, and K is a 
shape-dependent factor to account for non-uniform shear distribution in 
tte cross-section. The particular case of a rectangular cross-section 

3 

connesponds to I * bt /12 and A = bt, where t is the beam thickness and b 
its »idth. In order to simplify the algebra, it is convenient to combine 
th* lending and shear stiffness terms to obtain 



vhene a is the ratio EI/(kGA) with dimension length squared. For a 
rectangular cross-section, assuming K = 1 and incompressible material with 

2 

E * 3G, this ratio becomes a * t /4. 
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u (n+1) » u (n) 4 du (n+1) (4b) 

vhere the symbol ' is used to denote the perturbed quantities. The 
consistency of the algorithm is provided for in the computation of the 
right-hand-side of (4a), since the process is equivalent to the 
minimization of the residual 


» 


R u 


(n) 


which will be attained if u^ n \ * u, the exact value of the perturbed 

response. Stability is achieved if each displacement correction du^ 

is smaller (in an appropriate norm) than the preceeding term, du^ n \ Both 
conditions must be satisfied for the iteration to converge to the exact 
solution. 


The stability conditions can best be discussed in terms of an 
amplification matrix, which is derived next. Consider the form of 
equation (4a) in two consecutive iterations 

R du (n+1) « f - R u (n) Iteration (n) 


K du^ n ^ * f - R u (n ^ Iteration (n-1) 

and subtract the second from the first to obtain 


R du (n+1 > - R du (n > - -R(u (n > - u*"^) 

R du (n+1) - R du (n) - R du (n) (5) 


Premultiplication by R~* on both sides yields 


du (n+1) = (I - R " 1 R) du <n) 
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where 

i 

A = (I - K -1 K) (6) 

is the desired amplification matrix. The iteration will be stable if the 

dx/ n ^ decrease monotonically, i.e., if the spectral radius of every 
eigenvalue of A is less than unity. 

In this form, the determination of the spectrum of the amplification 
matrix A would require a considerable amount of computation. In general, 

the eigenvectors of the perturbed stiffness matrix K will be different 
from the eigenvectors of the unperturbed stiffness K. This will result in 
a nonsymmetric amplification matrix A. In addition, the size and 
structure of the amplification matrix in this form is entirely problem- 
dependent, so that it does not easily lend itself to analysis for the 
general case. 

Stability Analysis 

In order to circumvent some of the problems raised in the preceeding 
section, a Von Neumann stability analysis is performed on the difference 
pattern corresponding to the assembled system of equations at a typical 
internal node. Similar techniques have been used in studies of the 
stability of transient time integration schemes and nonlinear solution 
algorithms [1, 9]. 

The fundamental concept underlying these techniques is 
straightforward, even though the detailed derivations often require 
extensive algebraic manipulations. First, a set of stiffness equations 
corresponding to a typical node is extracted from the assembled stiffness 
equations. For a one-dimensional uniform mesh of tvo-noded beam elements, 
this will be a set of two equations, relating the shear and moment at node 
k to the translations and rotations at nodes k-1, k and k+1. Considering 
the linearly interpolated Timoshenko beam element in (3), these equations 


become 
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where and are the transverse shear force and moment at node k, and 
u^ and 0^ are the transverse displacement and the rotation at that node. 

In order to capture the characteristics of the assembled system for 
an arbitrary displacement vector, a sinusoidal displacement pattern of the 
form 



is imposed on the nodes of the one-dimensional mesh. Here, u and © are 
complex constants representing the relative magnitude and phase of the 
displacements, and u = 2nh/l where 1 is the (arbitrary) wavelength of the 
prescribed sinusoidal displacement pattern. Hence, the value « = 0 
corresponds to the two rigid-body modes (one translation and one 
rotation), and the value w = Jl will result in two displacement patterns 
which alternate signs between consecutive nodes. Any compatible 
displacement configuration of the discrete system may therefore be 
obtained by appropriate combination of a number of these "basic modes" 
with different u between 0 and n. 

Substituting (8) into (7) and using a few trigonometric identities, 
the following expression for the effective stiffness at an arbitrary « may 
be derived: 
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KGA 


r (1 - cosw) 
n 

-i sinw 


i sinw 

h /. x 2a 

j (1 + COSW) + jj- 


m 


» 



U 


i 

h 

(1 -cosw) 

4 


9 

fr 4 


(9) 


■ k(«) u 


This relation may be regarded as a "condensed" counterpart of the global 
stiffness equations, corresponding to a known displacement pattern (mode). 
Since no assumptions have been made on the value of w, the equation above 
must hold for all values of w that are compatible with the prescribed 
boundary conditions. 


The techniques outlined in the preceeding paragraphs may be used to 
construct a "condensed" counterpart of the algorithmic relation in (5) 
corresponding to a given value of w: 


k( w) du (n+1> = k(w) du (n) - k(w) du (n) (10) 

An amplification matrix relating consecutive displacement corrections for 
a given mode may be obtained by premultiplying (10) by k~*(w) to obtain 


du 


(n+1) 


a du 


(n) 


where 


a(w) = (I - k _1 (w) k(w)) (11) 

is the desired amplification matrix. Stability conditions associated with 
particular classes of perturbations of the one-dimensional beam mesh 
problem may then be derived from the study of the eigenvalues of (11). 

An interesting class of perturbation problems involves the (not 
necessarily uniform) elongation of the mesh. In the thin limit, the 
transverse shear constraint will impose the condition that 
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on the displacement solution of the unperturbed system. Hovever, in the 
perturbed (elongated) beam, a different transverse shear constraint is in 
effect, vhich is not satisfied by the displacement solution for the 
unperturbed system, i.e., 


dv 

dx 


* e 


From the form of (1), it is clear that a very large amount of shear 
deformation energy is generated when the displacement solution for the 
unperturbed problem is imposed on the perturbed system, even for seemingly 
"small" elongations of the mesh. 


In order to obtain a stability limit for this class of problems, a 
uniform elongation of the mesh is considered. The element length on the 

perturbed mesh thus becomes h = h(l+e), so that each element in the mesh 

is elongated by the same amount. It follows that k~*(w) may be expressed 
as 


u- 1 / h 

k («) - j; 


1 


(1 - cosw) 4 


h 2a 

2 (1+cosw) + jj- (1-cosw) - i sinw 


i sinw 


(1-cosw) 


( 12 ) 


and k(w) as 


k(w) 


RT7T) < 1 - cos “) 


-i sinw 


i sinw 

( 1+cosw) + (1-cosw) 


(13) 




The resulting amplification matrix will be 
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»(») 


2 2 

e c ,h . . „.h . . e ,h x 

l+c + 1+e < 4« ) cot ( 2 ) ic( S^ ) cot ( 2 ) _i 1+c < 2 ) cot( 2 ) 

2 

. g /h x .,«> e .h . ,2 y «x 

1 T7e ( Ii> cot( J> 171 - *W cot <5> 


and has eigenvalues of the fora 


(14) 


V“> ■ if E - 5 ITT (1 - J 


1 * ??» 


(15a) 


2 2 ^ 
x , x e 1 e f ,, , 4 

¥ w) = Ui ' 2 UT (1 + J 1 + ^ 


(15b) 


where 


.h . ..2.M. 

<4i> cot < 2 ) 


(15c) 


All values of « which are relevant to the analysis of the discrete 
system lie between 0 and n. The highest deformation modes representable 
by the discretized system correspond to » = it, which will result in 



Thus, even for relatively large t, the spectral radius of the 

amplification matrix will be less than unity and the corresponding 
deformation modes remain stable. At the opposite end of the spectrum lie 
the rigid body modes, corresponding to w = 0. As w approaches 0, the 

value of cot(^) becomes unbounded. This means that the rigid body modes 

are unconditionally unstable for any nonzero value of c, as expected. 

An asymptotic analysis of the eigenvalues of the amplification matrix 
for large values of c will yield 
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X 1 Sl 


*2 


c 2 f 2 

1+e 2 


From the behavior of cot(^), it can be concluded that f may become quite 
large for small values of ». This will result in | I > 1 even for 

seemingly small c, and will cause the associated deformation mode to grow 
unbounded. It should be emphasized, however, that the asymptotic 
expressions above typically represent reasonable approximations to the 
eigenvalues of i(u) only for values of c well above the stability limit 
and cannot be used to approximate the eigenvalues within the stability 
bounds . 

From the above discussion, it is clear that the higher deformation 
modes (with small values of w) will govern the stability of the algorithm. 
This is in contrast with the well known results for the stability of 
explicit time integration algorithms in dynamics, which are governed by 
the highest frequency modes present in the discretized system. Any 
attempts at, enhancing the stability of the perturbation algorithm must 
therefore take into account the fact that the displacement modes which 
require stabilization are among the most needed to represent the response 
of the perturbed system. Stability in the higher deformation modes must 
not compromise the accuracy of these modes, which rules out the use of 
conventional stabilization procedures. 

A Numerical Example 

A test problem was set up using a one-dimensional mesh of ten 
Timoshenko beam elements (NESSUS Element 98) with h * 2.00 and t * 0.25, 
and made of incompressible material (\> * 0.50). Three different cases 
were analyzed, corresponding to the following boundary conditions: 
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1. Cantilever beam. In this case, the lowest displacement mode has 

a wavelength of four times the beam length, corresponding to 

t^n/20. i 

2. Beam with both ends pinned. The wavelength of the lowest mode 
is twice the beam length, corresponding to c^n/10. 

3. Beam with both ends fixed. The wavelength of the lowest mode is 
equal to the beam length, corresponding to «^lt/5. 

The variation of the spectral radius of the eigenvalues of the 
amplification matrix in the lowest mode as a function of the perturbation 
parameter c is shown in Figures 3 to 5. In all cases, loss of stability 
was observed at the value of e corresponding to a spectral radius of 1.00, 
as predicted by the analysis. Similar behavior has been observed using a 
shell model (NESSUS Element 75) of the same problem. 

Conclusions 


The stability conditions for the elastostatic perturbation algorithm 
proposed in [2] have been described. A Von Neumann stability analysis of 
the algorithm was performed for the case of a uniform one-dimensional mesh 
of linearly interpolated Timoshenko beam elements. Closed form 
expressions for the stability limit in terms of the perturbation parameter 
e were derived for the case of uniform elongation of the mesh by a factor 
of (1+e). This form of perturbation has been observed to result in loss 
of stability with the current implementation of the algorithm even for 
seemingly small values of the perturbation parameter e. The stability 
limits predicted by the analysis are in full agreement with those observed 
in numerical experiments on one-dimensional meshes of linearly 
interpolated beam and plate elements. 

The development of general closed-form results for unstructured, 
multi-dimensional meshes subjected to non-uniform distortion does not 
appear to be practical. However, limited experience has indicated that 
the results for the uniform mesh can be used to obtain a conservative 
estimate to the stability limit for a more general mesh. Therefore, the 
development of "smart" algorithms to adaptively adjust the size of the 
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perturbation parameter in order to ensure convergence appears very 
promising. 
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Element 
h = 2.0 





Overall Length = 20.0 


FIGURE 1: Hesh for the One-Dimensional Model Problem 
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FIGURE 2: Typical Displacement Patterns for the Values of 

U Used in the Numerical Example 


Spectral Radius 
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Epsilon 


FIGURE 3: Spectral Radius of the Amplification Matrix as a 

Function of the Perturbation Parameter 
e for h/t = 8.00 and w = n/20 



Spectral Radius 
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Epsilon 


FIGURE A: Spectral Radius of the Amplification Matrix as a 

Function of the Perturbation Parameter 
e for h/t = 8.00 and w = n/10 


Spectral Radius 


Unstable 



Epsilon 


FIGURE 5: Spectral Radius of the Amplification Matrix as a 

Function of the Perturbation Parameter 
t for h/t = 8.00 and « = Jl/5 
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l.o DISCUSSION 

In private correspondence of January 13, 1986 with Y. N. Chen of the 
American Bureau of Shipping, we were informed that certain numerical algo- 
rithms in the Wu/FPI code lacked the precision to ensure "small" errors 
in resulting point probability estimates. Those subroutines for which 
improvements were suggested were: 

1. The normal CDF 

2. The inverse normal CDF 

3. The gamma function 

4. The shape parameter of the gamma function 

5. The EVD parameters 

ABS implimented improved numerical procedures in FPI and studied several 
examples. Because the differences in point probability estimates observed by 
ABS in old FPI and their new version seemed significant, a study was undertaken 

to carefully examine the approximate forms and to introduce improvements where 
appropriate. The improvement in the Euler constant (for EVD parameters) 

for 8 digit accuracy was trivial and was implimented immediately. Numerical 
algorithms for the other terms cited above were developed. Their performance 
was carefully examined. A detailed description of the approximation forms 

and their behavior is presented in Chapter 2. 

The forms presented in Chapter 2 were introduced into FPI, replacing 
their less accurate counterparts. FPI analysis using the old and new code 
was performed on several examples. The results are summarized in Chapter 3. 

Differences in the results of the old and new code are far less than 
observed by ABS in their version of the code. At this time, there is no 
explanation for the discrepancy 


2.0 APPROXIMATE FORMS OF FUNCTIONS USED IN PROBABILITY CALCULATIONS 

2.1 Gaasna Function 

Ref: Abramowitz, Handbook of Mathematical Functions. NBS. 


r(x) 



x > 0 


o 


The Asymptotic Formula 


£n T(x) - (x - j) £n x - 


* + 7 *” < 2 "> + ifc 


300x' 


1260x 5 1680x 7 


(x +■ • in |Atg x| < ir) 


After testing this formula, we found that when x — 6, ten digit accuracy is 
provided. If x is increased, this form is even more accurate. 

In this program, x is divided into two parts X - 6, and 0 < x < 6. 

If x- 6, use the asymptotic formula directly. If 0 < <6, then let 

N - INTEGER (x) 

Z - 6 - N + x 

and calculate £n T(Z) using the asymptotic formula. 


6-N 

Then let, £n T(x) - Zn T( Z) - [ £n(x + J - 1.0) 

J-l 


Example x - 1.9 


Z - 6. - 1 + 1.9 - 6.9 
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If more accuracy is needed, then increase the 6 in the above algorithm to 
7, 8, or a larger number 





Table 1. Performance of Gamma Function Approximation 
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Subroutine for Ganna Function 


DOUBLE PRECISION FUNCTION GAMMA(Y1,PI) 

IMPLICIT DOUBLE PRECISION (A-H,0-Z> 

X-YI+1.D+0 

Z*X 

I F ( X . GE . 6 . QD+E ) GO TO 456 
N=INT(X> 

Z=(6.0D+0>-N+X 
456 Y=l.D+0/Z**2 

ALG= < Z-. 5D+0) *DLOG <Z ) +. 5D+0*DLOG (PI *2. D+0) - 

* Z— ( 1 . D+0/ (12. D+0*Z ) )*( ( (Y/0. 14D+3-1 . D+0/0. 135D+3WY+ 

* 1 . D+0/. 3D+2> *Y— 1 . D+0) 

IF <X.GE.6.D+0)GO TO 457 
ITE=6-N 
DO 3 J=1 , ITE 
A=X+J-1 . D+0 
ALG=ALG-DLOG ( A ) 

CONTINUE 

7 GAMMA=DEXP (ALG) 

RETURN 
END 
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2.2 Biseccion Method for the Shape Parameter a of the Welbull Distribution 
The coefficient of variation in terms of the shape parameter a of 
the Weibull distribution is given as 

L x / r(l + 2a) 

Given C x , it is required to compute a. 

Define 

F(a) - - (1 + C 2 ) r 2 (l + a) + r(l + 2a) 

1 08 

Approximate ■ (C x ) • Then calculate F(a^), and let 
<* 2 * Qt 1 + Calculate F(a 2 > and let F12 - FCa^) * F(a 2 ). 

If F12 - 0; we know that the root will be bracketed by a and a 2> 

Then use the general bisection method as described below. 

If F12 > 0, there are four possible cases. 



Case 3 Case 4 
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If the function looks like Case 1 or Case 3, then let a 1 ■ <* 2 ; 


« 2 * ®2 + °* 


If the function looks like Case 2 or Case 4, then let ■ o^; °1 " °1 ” 

Then calculate F12 until F12 < 0, at which time the bisection method can be used 
(A) General Bisection Method 


1. a. 


a l + a 2 


(*) 


F13 - F(a x ) * F(a 3 ) 


If F13 < 0, 
If F13 > 0, 


If \a 1 - a 2 | > 


a 2 * a 3 


°1 " a 3 


10 ' go to (*) and repeat. 
10~ 7 STOP; Let a - a. 


(B) Performance 

Consider the Rayleigh Distribution 


C x - .522723201 



Asymptotic formula 

Exact 

a 

2.00000014531220 

2.0 
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Program: Bisection Method for Weibull Shape Parameter 


IMPLICIT DOUBLE PRECISION (A-H,0-2) 

F ( X , CO V , P I > *- < 1 . DB+CO V**2 ) *GAMMA ( X , P I > **2+GAMMA <2. *X ,PI > 
P I =4 . D0*DATAN < 1 . D0 ) 

COV*. 522723201 
X1=COV**<1.03) 

7 F1=F<X1 ,COV,PI> 

I F < DABS ( F 1 > . LE . 1 . D-7 ) GO TO 1 

X2=X 1 +• 1D0 

F2=F<X2,C0V,PI> 

F 1 2=F 1 *F2 

IF <F12. LT. B. ) GO TO 20 
IF <F1. GT. 0. . AND.F2.GT.F1 ) Xi=Xl-. 1D0 
IF(F1.LT.0. . AND.F2.GT.F1) X1=X2 
IF(F1.GT.0. . AND.F1.GT.F2) X1=X2 
IF(F1.LT.0. . AND.F1.GT.F2) Xi=Xi-. 1D0 
GO TO 7 

20 CONTINUE 

2 X3*(X1+X2)*.5D0 

F13=F<X1 ,COV,PI) #F<X3 t C0V,PI) 

IF (F13. LT. 0. ) X2=X3 
IF (F13. GT. 0. ) X 1 =X3 
DX=DABS ( X 1-X2) 

IF (DX . GE. 1 . D— 7) GO TO 2 
1 ALPHA= 1 . DB/ X 1 

WRITE (*,#) ' ALPHA = ALPHA 

ZZ=. 95D0 
DO 1000 1=1,21 
2Z=ZZ+. 05D0 

WRITEC*,*) ZZ , GAMMA (ZZ-1 . D0,PI ) 

1000 CONTINUE 

STOP 
END 



2.3 CDF of Normal Distribution ( 

Ref: Abramowitz: Handbook of Mathematical Functions. NBS 


rx 


P - *(x) 


i -it 2 


/ 2tt 


e 2 dt 


1 - Z(x) (b^t + b 2 t 2 + b 3 c 3 + b 4 t 4 + b 5 t 5 ) if x - 0 

Z(x)(bjt + b 2 t 2 + b 3 t 3 + b^t 4 + b 5 t 5 ) if x < o 


Where, 


i 1 2 

Z(x) - —j L. e " 2 x 


t - 


1 + px 


p - 0.231621 
b l - 0.319381530 

b 2 - -.356563782 

b 3 - 1.781477947 

b 4 - -1.821255978 

b $ « 1.330274420 


This approximation is advertised to produce errors in P of less than 10~^. 
(See P'tfonnance on Table 2, p. 15.) When x > 0, It appears that this level of- 
accuracy is being realized. For the very small P values associated with x < 0, 
errors are somewhat larger. But the operational range for structural reliabilit 
analysis is -5 < x < 5, and at worst we are getting four place accuracy. 

It is important to note that we cannot verify the accuracy of column (4) 
in Table 2, p. 15. During this investigation, some anomolies were discovered 


in the Abramowitz table. 
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Program: Standard Normal CDF 


PROGRAM CDFPDF 

IMPLICIT DOUBLE PRECISION <A-H t O-Z> 

COMMON /TWO/ PI,SPI2 
PI=4.D0*DATANC1.D0> 

SPI2=1. / (DSORT (2. D0*F'I ) ) 

X=- 1 1 . D0 
DO 1 1=1,22 
X=X+1.D0 
PH I =CDFNOR ( X ) 

XPHI=XINV (PHI > 

WRITE (*,*) X ,PHI , XPHI 

1 CONTINUE 
STOP 
END 

DOUBLE PRECISION FUNCTION CDFNOR(Z) 

C THIS FUNCTION COMPUTES THE NORMAL CDF. 

IMPLICIT DOUBLE PRECISION <A-H,0-Z> 

COMMON /TWO/ PI,SPI2 

DATA A/0. 31938153D0/ ,B/-0. 356563782D0/ ,C/ 1 . 781477937D3/ 
+D/-1 . 821 255 9 7800/ , E/’ 1 . 330274429D0/ 

EZ=- ( Z**2) *. 5D0 
CDFNOR=0. 0D0 

IF <EZ. LE. —200. 0D0) GO TO 1 
ZX=SPI2*DEXP (EZ ) 

IF (DABS <Z> .GT.6.D0) GO TO 2 
T=1 . D8/ ( 1 . D0+ (0. 2316419D0*DABS (Z) > ) 

CDFNOR=ZX*T* ( A+T* (B+T* (C+T* <D+T*E ) ) ) ) 

GO TO 1 

2 Z2-1.D0/ (Z*Z) 

CDFNOR=ZX* < 1 . D3-Z2+ ( 1 . D0-3. D3*Z2* ( 1 . D0-5. D3*Z2; ) > /DAE 
1 IF ( Z . GT. 0 . 0D0 ) CDFN0R=1 . 3D0-CDFNOR 

RE TLRn 
END 


C- 12 


2 . 4 Bisection Method for the Inverse Normal CDF, x»$ ^(p) 

Ref: Abraoowitz, Handbook of Mathematical Functions, NBS 


First, the following method is used to obtain. an approximation to x 

C + C t + C, t 2 
t ^ 1 2 


x x - * -1 (P) 


where, 


1 + d 1 t + d 2 t 2 + d 3 t 3 


t - /-2 in P, 


0 < P - .5 


C - 2.515517 
o 

C x - .802853 
C 2 - .010328 
d x - 1.432788 
d 2 * .18926 
d 3 - .001308 


This approximation gives only four digit accuracy. 

Define F(x, P) ■ P - 4>(x) . 

1. Let x^ = 4> ^(P) using the "crude" approximation above. Then 
F(x, P) looks like 


F(x 
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2. Let F x - F(x r P)* (A) 


If F x > 0, 


■ x^ + .001 


If F x < 0, 


x^ ■ x^ - .001 


If F x - 0, 


STOP 


Then in the second iteration, let 
F 2 - F(x 2 , P) 

Calculate F12 * F(x^, P) * F(x 2> P) 


If F12 < 0, use general Bisection Method 
If F12 > 0, then X l - x £ 


go back to (A) and repeat. 


The function $(•) is obtained using the form of Sec. 2.3. 



C— 14 



where F(x, P) * P - 4>(x) 

Fig. 2 Flow Chart for Inverse Normal Appi 



oximation 
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9.8659 E-10 



8.413447404 E-l 


9.772499379 E-l 


9.986500327 E-l 


9.999683139 E-l 


9.999999990 E-l 


- 2 . 00000000 


-.99999999 


3.61190816 E-ll 


.99999999 


2.00000000 


3.00000000 


3.99999999 


4.99999998 


7.9911351772922 


3.1671 E-5 


1.349898032 E-3 


2.275013195 E-2 


1.586552540 E-l 


5.0000 E-l 


8.413447460 E-l 


9.772498680 E-l 


9.986501019 E-l 


9.999683288 E-l 


Approximate form as 
described in Sec. 2.3 


The inverse is obtained 
using the column (2) 
values with the algorithm 
described in Sec. 2.4. 
Columns (1) and (3) 
should compare. 


Column (4) and 
column (1) should 
compare. See 
comments on p. 10, 










































Program: 


Bisection Method for Inverse Standard Normal CDF 


DCLBLE PRECISION FUNCTION XINV <Z) 

I -i- LICIT DOUBLE PRECISION (A-H.O-Z) 
F < X t P 1 ) =P 1 — CDFNOR < X ) 

Y=Z 

IF<Z.GT.0.5D0> Y= 1 . D0-Z 
IF (Z. EQ. 1 . D0) STOP 
C0*2.515517D0 
C1=0.B02853D0 
C2=0. 010328D0 
D1 =1 . 432788D0 
D2=0. 189269D0 
D3=0 . 00 1 308D0 

T= (—2. D0*DLOG ( Y) >**.5D0 
DNUM=C0+T* (C1+T*C2> 

DNOM= 1 . 0D0+T# (Dl+T* ( D2+T#D3 ) ) 

X=T- (DNUM/DNOM) 

IF < Z . LT .0. SD0) X=— X 
X1=X 

Fl=F(X t Z) 

80 IF <F1 . GT. 0. D0) X2=X1+.001D0 

IF (FI . LT .0. D0) X2=X 1— . 001D0 
IF (FI. EQ . 0 . D0 ) GO TO 2 
F2=F (X2,Z) 

FI I=F1*F2 

IF (F 12. LE.0.D0) GO TO 8 
X 1 =X2 
F 1 =F2 
GO TO 80 

3 X3= (X1+X2) *. 5D0 

F 1 3=F ( X 1 , Z ) *F ( X3 , Z > 

I F ( F 1 3 . LE . 0 . D0 ) X2=X3 
IF (FI 3. GT .0. D0) X 1 =X3 
DX=DABS < X 1-X2) 

IF (DX.GT. l.D-10) GO TO 8 
2 X INV=X 1 

RETURN 

end 
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3.0 EXAMPLES COMPARING OLD AND NEW FPI 

Following are several examples for which comparisons of results from 
old and new FPI are presented. These examples were those studied in an 
AME master's report by Jack T. L. Chang entitled, "Investigation of the Wu 
Algorithm for Computing Structural Reliability" (October 1985). In summary, 
introduction of the improved algorithms did not significantly alter the 
results, at least for the examples considered. 

Results for the improved ABS FPI program for those examples considered 
are given in parentheses. In Examples 4 and 5, the results of the Improved 
versions of the ABS and the UA codes differ significantly. There is at this 
time no explanation for the disagreement. An efficient Monte Carlo code 
for point probability estimates is under development. It will be able to 
check FPI calculations, but because the same numerical algorithms will be 
in both UA codes, the comparisons may not resolve this issue. 


EXAMPLE 1 


Note: ABS program results in parentheses 


-DATA 

FAILURE FUNCTION : g-R-(L + D) 


FAILURE EVENT : 9 < 0 


VARIABLE 

DISTRIBUTION 

MEAN /MEDIAN* 

STD. DEV. 

C.O.V. 

R 

WEIBULL 

50. 

5.0 

0.1 

L 

EVD 

10. 

2.0 

0.2 

D 

LOGNORMAL 

20. * 

3.034 

0.15 



ORIGINAL 

PR0GRAM (1) 

NEW 

program (2) 3 4 

DIFFERENCE 

R-F 

6 

2.768 

2.783 

(2.783) 

0.54 

Pf 

2.821 E-3 

2.6931 E-3 
(2.6931 E-3 

4.79 

_ 

Wu/FPI 

6 

2.692 

2.707 

(2.680) 

0.55 

Pf 

3.554 E-3 

3.398 E-3 
(3.682 E-3) 

i 

4.59 

Monte 
Carlo ^ 

Pf 

3.600 E-3 




(1) Developed by C. Kelly, Y. T. Wu 

(2) With improvements to numerical algorithms: (a) gamma function, (b) Weibull 

shape parameter, (c) standard normal cdf, (d) inverse normal cdf, 

(e) EVD parameters 

(3) Assumes new program is exact 

(4) Does not have improvements in numerical algorithms 
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EXAMPLE 2 Note: ABS program results in parentheses 


-DATA 

FAILURE Function : 



( 2N ) b + c’ 


( 2N ) c 


e 


s 


FAILURE EVENT : 9 < 0 


VARIABLE 

DISTRIBUTION 

MEAN /MEDIAN* 

STD. DEV. 

C.O.V. 

Cs 

EVD 

0.0015 

0.00015 

0.1 

. '1 

LOGNORMAL 

310.0 * 

145.10 

0.43 


LOGNORMAL 

9.14 * 

0.458 

0.05 



ORIGINAL 

PR0GRAM (1) 

NEW 

PROGRAM (2) 

DIFFERENCE 

Z (3) 4 

R-F 

6 

2.881 

2.881 

(2.881) 

0.00 

Pf 

1.981 E-3 

1.983 E-3 
(1.983 E-3) 

0.10 

Wu/FPI 

8 

2.851 

2.850 

(2.846) 

0.04 

Pf 

2.183 E-3 

2.183 E-3 
(2.215 E-3) 

0.00 

Monte 

Carlo 

Pf 

2.123 E-3 




(1) Developed by C. Kelly, Y. T. Wu 

(2) With improvements to numerical algorithms: (a) gamma function, (b) Weibull 

shape parameter, (c) standard normal cdf, (d) inverse normal cdf , 

(e) EVD parameters 

(3) Assumes new program is exact 

(4) Does not have improvements in numerical algorithms 





EXAMPLE 3 


Note : 


ABS program results In parentheses 


-DATA 

FAILURE FUNCTION : g - Xj + 2X 2 + 2X 3 +X A - 5( X 5 + X fi ) 
FAILURE EVENT : 9 < 0 


VARIABLE 

DISTRIBUTION 

MEAN /MEDIAN* 

STD. DEV. 

c.o.v. 

Xi 

LOGNORMAL 

119.4 

* 

12. 

0.1 

X 2 

LOGNORMAL 

119.4 

* 

12. 

0.1 

X3 

LOGNORMAL 

119.4 

* 

12. 

0.1 

Xh 

LOGNORMAL 

119.4 

* 

12. 

0.1 

X 5 

LOGNORMAL 

38.31 

* 

12. 

0.3 

XS 

LOGNORMAL 

47.89 

* 

15. 

0.3 



ORIGINAL 

PROGRAM (1) 

NEW 

PR0GRAM (2) 3 

DIFFERENCE 

R-F 

8 

2.348 

2.348 

(2.348) 

0.00 

Pf 

0.942 E-2 

0.943 E-2 
(0.943 E-2) 

0.1 

Wu/FPI 

e 

2.235 

2.234 

(2.256) 

0.04 

pf 

1.274 E-2 

1.274 E-2 
(1.204 E-2) 

0.00 

Monte 
Carlo (4) 

1 

pf 

1.221 E-2 




(1) Developed by C. Kelly, Y. T. Wu 

(2) With improvements to numerical algorithms: (a) gamma function, (b) Weibull 

shape parameter, (c) standard normal cdf, (d) inverse normal cdf , 

(e) EVD parameters 

(3) Assumes new program is exact 

(4) Does not have improvements in numerical algorithms 





C-21 


EXAMPLE 4 Note: ABS program results in parentheses 


-DATA 

FAILURE FUNCTION : g - X + X + X + X + X - Y - Y - Y - Y - Y 

1234512345 

FAILURE EVENT : 9 < 0 


VARIABLE 

DISTRIBUTION 

MEAN /MEDIAN* 

STD. DEV. 

c.o.v. 

Xl 

WEIBULL 

10.0 

3.5 

0.35 

X2 

WEIBULL 

10.0 

3.5 

0.35 

x 3 

WEIBULL 

10.0 

3.5 

0.35 

X4 

WEIBULL 

10.0 

3.5 

0.35 

x 5 

WEIBULL 

10.0 

3.5 

0.35 

Yl 

EVD 

5.0 

1.75 

0.35 

Y2 

EVD 

5.0 

1.75 

0.35 

y 3 

EVD 

5.0 

1.75 

0.35 

y 4 

EVD 

5.0 

1.75 

0.35 

YS 

EVD 

5.0 

1.75 

0.35 



ORIGINAL 

PROGRAM (1) 

NEW 

program (2) 3 4 

DIFFERENCE 

R-F 

8 

2.945 

2.959 

(2.959) 

0.47 

Pf 

1.615 E-3 

1.545 E-3 
(1.545 E-3) 

4.53 

Wu/FPI 

8 

2.866 

2.877 

(2.810) 

0.38 

p f 

2.078 E-3 

2.011 E-3 
(2.477 E-3) 

3.33 

Monte 
Carlo <4) 

Pf 

2.140 E-3 




(1) Developed by C. Kelly, Y. T. Wu 

(2) With improvements to numerical algorithms: (a) gamma function, (b) Weibull 

shape parameter, (c) standard normal cdf, (d) inverse normal cdf , 

(e) EVD parameters 

(3) Assumes new program is exact 

(4) Does not have improvements in numerical algorithms 





EXAMPLE 5 


Note: ABS program results in parentheses 


•DATA 

FAILURE FUNCTION ! 9 - X, + X + X + X + X - Y 

1 2 3 4 5 1 

FAILURE EVENT : 9 < 0 



VARIABLE 

DISTRIBUTION 

ME AN /MEDIAN* 

STD. DEV. 

C.O.V. 

Xl 

EVD 

10.0 

4.0 

0. 4 

X 2 

WEIBULL 

10.0 

4.0 

0* 4 

x 3 

LOGNORMAL 

9.2847 * 

4.0 

0. 4 

Xu 

EVD 

10.0 

4.0 

0. 4 

Xs 

WEIBULL 

10.0 

4.0 

0. 4 

Y 1 

EVD 

5.0 

2.0 

0.4 

Y 2 

WEIBULL 

5.0 

2.0 

0. 4 i 

X 3 

LOGNORMAL 

4.6424 * 

2.0 

0. 4 


EVD 

5.0 

2.0 

0.4 

Y5 

WEIBULL 

5.0 

2^0 

0.4 



ORIGINAL 

PR0GRAM (1) 

NEW 

program (2) 3 4 

DIFFERENCE 

z (3) 

R-F 

B 

2.649 

2.652 
(2. 652) 

0.11 

Pf 

4.031 E-3 

4.003 E-3 
(4.003 E-3) 

0.70 

Wu/FPI 

B 

2.696 

2.698 

(2.658) 

0.07 

Pf 

3.508 E-3 

3.491 E-3 
(3.984 E-3) 

■ 

0.49 

Monte 
Carlo < 4 > 

Pf 

3.643 E-3 




(1) Developed bv C. Kelly, Y. T. Wu 

(2) With improvements to numerical algorithms: (a) gamma function, (b) Weibull • 

shape parameter, (c) standard normal cdf, (d) inverse normal cdf , 

(e) EVD parameters 

(3) Assumes new program is exact 

(4) Does not have improvements in numerical algorithms 





C-23 


EXAMPLE 6 Note: ABS program results in parentheses 


FAILURE FUNCTION : 9 - R - A-) 

inr 

FAILURE EVENT : 9 < 0 


VARIABLE 

DISTRIBUTION 

MEAN /MEDIAN* 

STD. DEV. 

C.O.V. 

R 

NORMAL 

170. 

25. 

0.14706 

D 

NORMAL 

29.4 

3 * 

0.10204 



ORIGINAL 

PROGRAM 

NEW 

program 1 (2) 3 

DIFFERENCE 

R-F 

6 

2.902 

2.902 

(2.902) 

0.00 

Pf 

1.856 E-3 

1.856 E-3 
(1.856 E-3) 

0.00 

Wu/FPI 

6 

2.835 

2.834 

(2.833) 

0.03 

Pf 

2.296 E-3 

2.297 E-3 
(2.306 E-3) 

0.04 

Exact 

solution 

1 

Pf 

2.301 E-3 




(1) Developed by C. Kelly, Y. T. Wu 

(2) With improvements to numerical algorithms: (a) gamma function, (b) .Weibull 

shape parameter, (c) standard normal cdf, (d) inverse normal cdf , 

(e) EVD parameters 

(3) Assumes new program is exact 





































C- 25 


EXAMPLE 7 Note: ABS program results In parentheses 


•DATA 

FAILURE FUNCTION : 9 - R - / 3 qop 2 + 1.92T 2 

FAILURE EVENT : 9 < 0 


VARIABLE 

DISTRIBUTION 

MEAN/MEDIAN* 

STD. DEV. 

C.O.V. 

R 

WEIBULL 

48.0 

3.0 

0.0625 

P 

LOGNORMAL 

0.987 * 

0.16 

0.16 

T 

EVD 

20.0 

2.0 

0.1 



ORIGINAL 

PROGRAM^ 1 * 

NEW 

PROGRAM* 2 * 

DIFFERENCE 

Z* 3 4 > 

R-F 

6 

3.094 

3.085 
f 3 . 0 8 5 

0.29 

Pf 

0.988 E-3 

1.018 E-3 
(1.016 E-3) 

2.95 

Wu/FPI 

6 

2.893 

2.886 

(2.868) 

0.24 

Pf 

1.911 E-3 

1.950 E-3 
(2.064 E-3) 

2.0 

Monte 

C.rlo <4> 

Pf 

1.800 E-3 




(1) Developed by C. Kelly, Y. T. Wu 

(2) With improvements to numerical algorithms: (a) gamma function, (b) Weibull 

shape parameter, (c) standard normal cdf, (d) inverse normal cdf, 

(e) EVD parameters 

(3) Assumes nev program is exact 

(4) Does not have improvements in numerical algorithms 





EXAMPLE 8 


-DATA 


FAILURE FUNCTION 
FAILURE EVENT 


DISTRIBUTION 

MEAN /MEDIAN* 

WEIBULL 

EVD 

LOGNORMAL 

150. 

100. 

0.1 * 


STD. DEV. 


0.1414 


0.16667 

0.2 

1.0 


ORIGINAL NEW j DIFFERENCE 

PROGRAM^ PROGRAM^ 


R-F 

8 

2.060 

2.067 

0.34 

Pf 

1.968 E-2 

1.938 E-2 

1.55 

Ur/FPT 

8 

1.967 

i 

i 

1.974 

0.35 

n u / r L 1. 

Pf 

2.461 E-2 

2.419 E-2 

1.74 

Monce 

Carlo 

Pf 

2.412 E-2 




(1) Developed by C. Kelly, Y. T. Wu 

(2) With improvements to numerical algorithms: (a) gamma function, (b) Weibull 

shape parameter, (c) standard normal cdf, (d) inverse normal cdf , 

(e) EVD parameters 

(3) Assumes new program is exact 

(4) Does not have improvements in numerical algorithms 

































EXAMPLE 8 


-DATA 

FAILURE FUNCTION : 9 - K - S /wA 
FAILURE EVENT : g < 0 


VARIABLE 

DISTRIBUTION 

MEAN /MEDIAN* 

STD. DEV. 

C.O.V. 

K 

WE I BULL 

150. 

25.0 

0.16667 

S 

EVD 

30. 

16.0 

0.2 

A 

LOGNORMAL 

0.1 * 

0.1414 

1.0 


u 

s 

a 

s 



ORIGINAL 

PR0GRAM (1) 2 3 4 

NEW 

PROGRAM 

DIFFERENCE 

R-F 

s 

2.482 

2i. 490 

0.32 

Pf 

6.534 E-3 

6.382 E-3 

2.38 

Wu/FPI 

6 

2.380 

2.389 

0.38 

Pf 

8.672 E-3 

8.453 E-3 

2.59 

Monte 
Carlo < 4 > 

Pf 

8.630 E-3 




(1) Developed bv C. Kelly, Y. T. Wu 

(2) With improvements to numerical algorithms: (a) gamma function, (b) Veibull 

shape parameter, (c) standard normal cdf, (d) inverse normal cdf , 

(e) EVD parameters 

(3) Assumes new program is exact 

(4) Does not have improvements in numerical algorithms 










C-28 


EXAMPLE 8 Note: ABS program results In parentheses 


-DATA 

FAILURE FUNCTION : g - K - S /Ha 
FAILURE EVENT : g < 0 


VARIABLE 

DISTRIBUTION 



K 

S 

A 

WE I BULL 
EVD 

LOGNORMAL 

otAii / nt u i an K 

150. 

60- 

0.1 * 

2>ID. DEV. 

25.0 

12.0 
0.1414 

C.O.V. 

0.16667 

0.2 

1.0 




ORIGINAL 

PROGRAM^ 

NEW 

PROGRAM 1 (2) 

DIFFERENCE 

* (3) 4 

R-F 

8 

3.006 

3.018 

£3.018) 

0.40 


Pf 

1.323 E-3 

1.272 E-3 
(1.272 E-ll 

4.01 

Wu/FPI 

8 

2.892 

2.905 
(2. 897) 

0.45 

Pf 

1.914 E-3 

1.835 E-3 

(1.886 E-1) 

4.30 

Monte 
Carlo (4 > 

Pf 

1.870 E-3 



— 






(1) Developed bv C. Kelly, Y. T. Wu 

(2) Wthl.provM.nt. to nuMrlcU algorithms: (a) gamma function, (b) Wa-bull 

shape paramatar, (c) standard normal cdf, (d) Invars, normil cdf 
(e) EVD parameters ’ 

(3) Assumes new program Is exact 

(4) Does .not. have improvements in numerical algorithms 
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EXAMPLE 9 Noce: ABS program results in parentheses 


-DATA 

FAILURE FUNCTION : 9 - A - N 0 


fop 

G <Y Ae Q ) Y 


+ 


1 ~ fpp 
H (Y Ae 0 ) n 


} 


FAILURE EVENT : 9 < 0 


VARIABLE 

DISTRIBUTION 

MEAN /MEDIAN* 

STD. DEV. 

C.O.V. 

A 

LOGNORMAL 

1.0 

* 

0.3132 

0.3 

fpp 

NORMAL 

0.7 


0.07 

0.1 

G 

LOGNORMAL 

0.222 

* 

0.0956 

0.4 

Y 

LOGNORMAL 

1.0 

* 

0.1517 

0.15 

Ae o 

EVD 

0.0005 


0.00008 

0.16 

H 

LOGNORMAL 

1.673 

* 

0.7208 

0.4 



ORIGINAL 

PR0GRAM (1) 

NEW 

program (2) 3 4 

DIFFERENCE 

Z< 3) 

R-F 

6 

2.384 

2.385 

(2.385) 

0.04 

Pf 

8.552 E-3 

8.550 E-3 
(8.550 E-3) 

0.02 

Wu/FPI 

S 

2.338 

2.338 

(2.315) 

0.00 

Pf 

9.696 E-3 

9.691 E-3 
(10.320 E-3) 

0.05 

Monte 

<W 4) 

Pf 

• 

i 

10.020 E-3 




(1) Developed by C. Kelly, Y. T. Wu 

(2) With improvements to numerical algorithms: (a) gamma function, (b) Weibull 

shape parameter, (c) standard normal cdf, (d) inverse normal cdf, 

(e) EVD parameters 

(3) Assumes new program is exact 

(4) Does not have improvements in numerical algorithms 





EXAMPLE 10 


-DATA 


FAILURE FUNCTION : 9 * R - Xj - X 2 
FAILURE EVENT : g<0 


VARIABLES 

Xi , 1 - 1 , 2 . 

DISTRIBUTION 

All XI are Chi-Square distribution 
with degree of freedom v ■ 1 . 

MEAN 

1.0 

STD. DEV. 

1.4142 

C. 0. V. 

1.4142 

CONSTANT , R 

3,4,5. 



ORIGINAL 

PR0GRAM (1) 2 3 

NEW 

program (2J 

DIFFERENCE 

Z< 3 > 

R-F 

6 

2.584 

2.583 

0.04 

Pf 

0.489 E-2 

0.490 E-2 

0.20 

Wu/FPI 

6 

2.178 

2.178 

0.00 

Pf 

1.471 E-2 

1.471 E-2 

0.00 

Exact 
solut ion 

Pf 

1 

1.110 E-2 




(1) Developed by C. Kelly, Y. T. Wu 

(2) With improvements to numerical algorithms: (a) gamma function, (b) Weibull 

shape parameter, (c) standard normal cdf, (d) inverse normal cdf, 

(e) EVD parameters 

(3) Assumes new program is exact 





C-31 



ORIGINAL 

PROGRAM (1) 

NEW 

PROGRAM ( 2 

DIFFERENCE 

* (3) 

R-F 

0 

3.676 

3.675 

0.03 

Pf 

1.186 E-4 

1.189 E-4 

0.08 

Wu/FPI 

6 

3.393 

3.390 

0.09 

Pf 

3.456 E-4 

i 

3.494 E-4 

1.09 

Monce 
Carlo (4) 

l 

Pf 

3.350 E-4 




\. 



ORIGINAL 

PROGRAM (1) 

NEW 

PROGRAM^ 

DIFFERENCE 

Z< 3 > 

R-F 

0 

4.735 

4.735 

0.00 

Pf 

1.096 E-6 

1.098 E-6 

0.18 

Wu/FPI 

0 

4.545 

4.535 

0.22 

Pf 

2.745 E-6 

2.879 E-6 

4.65 

Exact 

solution 

Pf 

3.730 E-6 


















































EXAMPLE 10 


-DATA 


FAILURE FUNCTION 


FAILURE EVENT 


VARIABLES 


DISTRIBUTION 


MEAN 


STD. DEV. 


C 


CONSTANT , R 


COMPARISONS OF SAFETY INDEX AND PROBABILITY OF FAILURE , Pf 


xi , 1 » 1, 2. 3, 4, 5. 


All Xi are Chi— Square distribution 
vith degree of freedom v - 1 . 


1.0 


1.4142 


1.4142 


3,4,5. 


ORIGINAL NEW 

PROGRAM^ PROGRAM^ 


2.049 


R-F 

B 

2.049 

Pf 

2.023 E-2 

Wu/FPI 

B 

1.302 

Pf 

' 9.652 E-2 

Exact 

solution 

Pf 

1.090 E- ; 


1.301 


DIFFERENCE 

*(3) 


0.08 


0.03 

































ORIGINAL 

PROGRAM* 

NEW 

PROGRAM* 2 ^ 

DIFFERENCE 

Z* 3 > 

R-F 

6 

3.241 

3.241 

0.00 

Pf 

5.954 E-4 

5.966 E-4 

0.20 

Wu/FPI 

6 

2.447 

2.447 

0.00 

Pf 

7.220 E-3 

7.195 E-3 

0.38 

Exact 

solution 

! 

pf 

6.840 E-3 




I . 



ORIGINAL 

PROGRAM* 

NEW 

PROGRAM* 2 ^ 

DIFFERENCE 

Z* 3) 

R-F 

6 

4.380 

4.380 

0.00 

Pf 

5.930 E-6 

5.951 E-6 

0.35 

Wu/FPI 

B 

3.574 

3.578 

0.11 

Pf 

1.761 E-4 

1.733 E-4 

1.62 

Exact 

solution 

Pf 

1.390 E-4 
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EXAMPLE 10 
-DATA 

FAILURE FUNCTION : g - R -X r X 2 -X 3 -X 4 -X 5 -X 6 -X 7 -X 8 -X 9 -X l0 


FAILURE EVENT : g < 0 


VARIABLES 

XI , i - 1,2,3,4,5,6,7,8,9,10. 

DISTRIBUTION 

All XI are Chi-Square distribution 
with degree of freedoa v ■ 1 . 

MEAN 

1.0 

STD. DEV. 

1.4142 

C. 0. V. 

1.4142 

mmm 

4,5,6. 


-COMPARISONS OF SAFETY INDEX AND PROBABILITY OF FAILURE , Pf 



ORIGINAL 

PROGRAM (1) 

NEW 

PROGRAM^ 

DIFFERENCE 

R-F 

6 

2.595 

2.595 

0.00 

Pf 

4.733 E-3 

4.725 E-3 

0.17 

Wu/FPI 

6 

1.254 

1.254 

0.00 

Pf 

1.049 E-l 

1.050 E-l 

0.09 

Exact 

solution 

Pf 

i 

0.060 E-2 
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EXAMPLE 10 



Exact 

solution 


Pf 


2.010 E-l 
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EXAMPLE 11 


-DATA _ 

FAILURE FUNCTION : g • 2.5 - N — H log * p 

1 + e o p o 

FAILURE EVENT : g < 0 


variable 

DISTRIBUTION 

MEAN/MEDIAN* 

STD. DEV. 

C.O.V. 

N 

NORMAL 

1.0 

0.1 

0.10 

Cc 

NORMAL 

0.396 

0.099 

0.25 

# 0 

NORMAL 

1.19 

0.1785 

0.15 

H 

NORMAL 

168.0 

8.4 

0.05 

Po 

NORMAL 

3.72 

0.186 

0.05 

AP 

NORMAL 

0.35 

0.07 

0.20 



ORIGINAL 

PROGRAM (1) 

NEW 

program (2) 3 4 

DIFFERENCE 

R-F 

6 

2.439 

2.439 

0.00 

Pf 

7.363 E-3 

" 363 E-3 

0.00 

Wu/FPI 

8 

2.499 

2.499 

0.00 

Pf 

6.235 E-3 

6.229 E-3 

0.10 

Monte 
Carlo (4 > 

Pf 

6.330 E-3 




(1) Developed by C. Kelly, Y. T. Wu 

(2) With improvements to numerical algorithms: (a) gamma function, (b) Weibull 

shape parameter, (c) standard normal cdf, (d) inverse normal cdf, 

(e) EVD parameters 

(3) Assumes new program is exact 

(4) Does not have improvements in numerical algorithms 





EXAMPLE 12 


-DATA 

FAILURE FUNCTION : g 


N C L B 


3/2 


FAILURE EVENT j g < 0 


RQl 


VARIABLE 


DISTRIBUTION 


MEAN /MEDIAN* 


STD. DEV. 


C.O.V. 


N 

C 

L 

H 

R 

Qt 


NORMAL 

NORMAL 

NORMAL 

NORMAL 

NORMAL 

EVD 


1.0 

3.85 

93.4 

15.0 

0.7 

9146.0 


0.2 

0.2695 

5.604 

0.9 

0.098 


3201.1 


0.20 

0.07 

0.06 

0.06 

0.14 

0.35 




ORIGINAL 

PROGRAM^ 15 

NEW 

PROGRAM^ 

DIFFERENCE 

R-F 

e 

2.715 

2.715 

0.00 

Pf 

3.309 E-3 

3.315 E-3 

0.18 

Wu/FPI 

6 

2.651 

2.651 

0.00 

Pf 

4.019 E-3 

1 

4.017 E-3 

0.05 

Monte 
Carlo 1 2 3 (4) 

Pf 

4.043 E-3 




(1) Developed by C. Kelly, Y. T. Wu 

(2) With improvements to numerical algorithms; (a) gamma function, (b) Weibull 

shape parameter, (c) standard normal cdf, (d) inverse normal cdf , 

(e) EVD parameters 

(3) Assumes new program is exact 

(4) Does not have improvements in numerical algorithms 
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EXAMPLE 13 


-DATA 

FAILURE FUNCTION : g ■ A + B H 3 + c Y S 3 + D Z Q 3 
FAILURE EVENT : g < 0 


VARIABLE 


DISTRIBUTION 


MEAN/MEDIAN* 


STD. DEV. 


C.O.V. 


X 

Y 

Z 

R 

S 

A. 


WE I BULL 
EVD 

LOGNORMAL 

EVD 

LOGNORMAL 
WE I BULL 


10.0 

5.0 

9.5782 * 

10.0 

4.7891 * 

10.0 


3.0 

0.30 

1.5 

0.30 

3.0 

0.30 

3.0 

0.30 

1.5 

0.30 

3.0 

0.30 




ORIGINAL 

PROGRAM (1) 2 3 

NEW 

PROGRAM^ 

DIFFERENCE 

R-F 

B 

2.625 

2.631 

0.23 


Pf 

4.327 E-3 

4.252 E-3 

1.76 

Wu/FPI 

S 

2.720 

2.724 

0.15 

Pf 

3.269 E-3 

3.223 E-3 

1.43 

Monte 
Carlo (4) 

Pf 

3.357 E-3 




(1) Developed by C. Kelly, Y. T. Wu 

(2) With improvements to numerical algorithms: (a) gamma function, (b) Weibull 

shape parameter, (c) standard normal cdf, (d) inverse normal cdf, 

(e) EVD parameters 

(3) Assumes new program is exact 

(4) Does not have improvements in numerical algorithms 





EXAMPLE 13 


FAILURE FUNCTION : g-A + BXR 1 2 3 +CYS 4 + DZQ ! 


FAILURE EVENT : g < 0 


VARIABLE 

DISTRIBUTION 

MEAN/MED IAN* 

STD. DEV. 

C.O.V. 

X 

WEIBULL 

10.0 

3.0 

0.30 

Y 

EVD 

3.0 

1.3 

0.30 

Z 

LOGNORMAL 

9.5782 * 

3.0 

0.30 

R 

EVD 

10.0 

3.0 

0.30 

S 

LOGNORMAL 

4.7891 * 

1.5 

0.30 

Q 

WEIBULL 

10.0 

3.0 

0.30 


a 

6 

Y 



ORIGINAL 

PROGRAM* 15 

NEW 

PROGRAM* 25 

DIFFERENCE 

Z* 35 

R-F 

6 

2.290 

2.290 

0.00 

Pf 

1.102 E-2 

1.102 E-2 

0.00 

Wu/FPI 

0 

2.410 

2.422 

0.50 

Pf 

7.983 E-3 

7.720 E-3 

3.41 

Monte 

Carlo* 45 

Pf 

8.020 E-3 




(1) Developed by C. Kelly, Y. T. Wu 

(2) With improvements to numerical algorithms: (a) gamma function, (b) Weibull 

shape parameter, (c) standard normal cdf, (d) inverse normal cdf , 

(e) EVD parameters 

(3) Assumes new program is exact 

(4) Does not have improvements in numerical algorithms 
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EXAMPLE 13 


-DATA 

FAILURE FUNCTION 


J • A + B n 5 + c Y ^ + DZQ 5 


FAILURE EVENT : g < 0 


VARIABLE 

DISTRIBUTION 

\ KEAN/MEDIAN* 

STD. DEV. 

C.O.V. 

X 

WEIBULL 

10.0 

3.0 

0.30 

Y 

EVD 

5.0 

1.5 

0.30 

Z 

LOGNORMAL 

9.5782 * 

3.0 

0.30 

R 

EVD 

10.0 

3.0 

0.30 

S 

LOGNORMAL 

4.7891 * 

1.5 

0.30 

0 

WEIBULL 

10.0 

3.0 

0.30 



ORIGINAL 

PR0GRAM (1) 

NEW 

procram (2) 3 4 

DIFFERENCE 

% (3) 

R-F 

6 

2.388 

2.392 

0.17 

Pf 

8.478 E-3 

8.369 E-3 

1.30 

1 

Wu/FPI 

6 

2.546 

2.549 

0.12 

Pf 

5.451 E-3 

5.405 E-3 

0.85 

Monce 

Carlo (4 > 

1 

Pf 

5.776 E-3 




(1) Developed bv C. Kelly, Y. T. Wu 

(2) With improvements to numerical algorithms: (a) gamma function, (b) Weibull 

shape parameter, (c) standard normal cdf, (d) inverse normal cdf , 

(e) EVD parameters 

(3) Assumes new program is exact 

(4) Does not have improvements in numerical algorithms 
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1.0 INTRODUCTION 

1.1 Some Definitions and Preliminary Remarks 

Let Y denote the response variable . Assume that Y will be a function 
of the random vector £ of design factors 

where Y - f(g) (1.1) 

X - (x x , x 2 , . . . X R ) 

This function is explicit and defined only through the data base generated 
by NESSUS, 

(Y t ; i - 1, J (1.2) 

where J is the number of solution points. 

It will be assumed that the basic statistical parameters for each X^ will be 
the mean and standard deviation , denoted as, 

E(X.) = p. V(X )’- o 2 ± (1.3) 

The vector parameter for is defined as, 

9 i - (u i , o i ) (1.4) 

And the parameter for X is a vector of K elements 0^, for a total of 2K statistical 
parameters 

e - (6., e 2 , . . . e K ) (1 * 5) 

As input to FPI , the statistical distributions of each X A must be specified. 


This includes the values of 0 . 


Consider a random sample of X ± ; j « l,n. The sample size is n. 


The estimators of y ^ and o^^ are, 


, n 

“i-n Ih.i 


( 1 . 6 ) 


n 


’i * ET * <X 1,J - “!>' 
j-1 


(1.7) 


The estimated parameters for X^ are 


e i ■ <*i. °i> 


( 1 . 8 ) 


and for all 


6 ‘ ( 6 1 , ® 2 , . . V 


(1.9) 


Using G, FPI constructs the distribution function (cdf) of Y, F^(y , G) . 
This is an estimate of the underlying cdf F^Cy), . . . the function which 
nature has chosen. An illustration of F^ is provided in Fig. 1.1. 

The distribution parameters G used to construct F^ are based on random 
samples. But the estimators Q are random variables themselves. There is 
uncertainty in the parameters which is reflected in F y . This uncertainty 
can be described by error bounds (or confidence intervals) as illustated 

in Fig. 1.1. It is the goal of this analysis to develop an operational pro- 
cedure for efficient estimation of these error bounds for implimentation 
in FPI. 

In classical statistics, 6 is considered to be chosen by nature and is 
a real number whose value remains forever unknown. 


The estimators 6 are 
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random variables and are used for point estimates of 6 and for contracting 
confidence intervals on 6. But for ease of analysis of confidence bounds, 
it is often convenient to use a "Bayesian approach" in which 6 is considered 
as a random variable reflecting the fact that its value is uncertain. Con- 
tinuing this role reversal the estimators 0 are assumed to be constant, 
and equal to the expected values of G. The value of this approach lies in the 
fact that if one can establish the distribution of 0, then upper and lower 
confidence bounds are just the appropriate percentage points. 

As an example, let the mean y of a normal variate be a random variable 
having a mean of y and standard deviation of a/Su*. A direct computation of 
the upper 95% and lower 5% points produces the 90% confidence interval on y. 
While deviating from classical statistics, this approach has experienced 
increased popularity in recent years. 
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1.2 Statement of the Problem 

Consider again the cdf of Y as shown in Fig. 1.2. If 0 were the actual 
values of 0 which nature has chosen, then we have perfect knowledge of the 
inherent variability of Y. would define precisely the distribution 
of Y. But if 0 is a random variable, then for a given F^, say F^, there will 
be uncertainty in the value of Y which produces F'. Thus, Y will be a random 
variable. 

It is important to note that what is really wanted is not the uncertainty 
of Y given F^, but rather that of F^ for a given Y, say y*. Thus, the general 
goal of this analysis will be to develop a practical algorithm for computing 
the error (or confidence) bounds on F^ for a given Y ■ y'. 

1.3 Response Variable as a Function of the Parameters 

To define the distribution of Y it is necessary to have an explicit 
expression for Y in terms of X. This function is constructed from the data 

i 

base (Eq. 1.2) as a polynomial. 

K K 

Y - f(X) - a +1 a X + l b xj + £ c X X (1.10) 

i«l 1 i-1 i,j 1 1 J 

i^J 

Now consider the distribution of X^, defined by the cdf, F^x; 

* 

and shown in Fig. 1.3. Let be the design point value corresponding to 

* 

y . FPI computes a design point JK when computing F^Cy’), and in fact, must 
satisfy, 

y ' * f(x*) (i.il) 

* * * 

The cdf corresponding to X^ is denoted as F^. At F^, X^ can be written in 
terms of by inverting the cdf 


W “ F i 1 (F i } 


( 1 . 12 ) 




X 


Fig. 1.3. The cdf of and the pdf of given F,. resulting from 
uncertainty in ©.. 


0-3 
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Upon substitution of each X^ into Eq. 1.10, Y can now be expressed in 
terms of 


Y - g(0) (1.13) 

1.4 Distribution of the Response Variable at a Given 

Because 0 is a random variable, Y is a random variable. And because 

the uncertainty of each was derived at F* (and X*), it follows that Eq. 1.10 

defines the distribution of Y at F'. Let Y denote the random variable, 

o 

Y at F'. The mean of Y should be "close to" y'. 

o 

Let the cdf and pdf of Y^ be denoted as F^ and f^ respectively; and let 

the mean and standard deviation of Y be u and a . 

o o o 


W i (0) - E(Y) - E[g(0) ] - y' (1.14) 

o*(0) - V(Y) - V[g(0) ] (1.15) 

1.5 Confidence Bounds on Y 
o 

Let a denote the confidence level, and let y T and y 1T denote the upper 

i. u 

and lower confidence bounds. These terms are related by the probability 
expression, 

P[y L ' Y o • y U J " ° 


And it follows that, 


p tv 0 - y L ) ■ f 0 (y L : y'. c 0 > ■ 1 -j- 2 


p tv o iy 0 ] - r< v y\ c o ) - — ° 


2 
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The upper and lover bounds are 



y and y define points on the confidence boundard as shown In Fig. 1.4. 

L U 

Translating horizontal confidence bound to a vertical bound statement, 

p [Y, > y D l • P[F < F’lyy]. 


Thus, one point on the lover confidence bound of at y^ is obtained. 
Similarly, 

P[Y o * y L ] - P[F > P ’l y L 1 

And a point on the upper confidence bound of F^ at y^ is defined. 

In general, then the confidence boundaries would have to be constructed 
on a point by point basis using several values of y f . 

A simpler scheme for/ estimating the error bounds for F^ at y’ using 
calculations at y' only will be presented in the next chapter. 
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2.0 EXAMPLE: "FIRST ORDER" ERROR BOUNDS 

2.1 Preliminary Remarks 

The problem of constructing error bounds on the cdf of response variable 
Y, as described in Sec. 1.0, may be too general to be practical. An example 
provided In this section Illustrates how an approximation to .the error bounds 
can be constructed using an algorithm which is simple enough to be included 
without great difficulty (we hope) in a probabilistic structural code.’ 

2.2 The Response Variable, Y 

Assume that Y is linear in £ in the neighborhood of y'. 

K 

Y - a + T a X (2.1) 

0 1-1 1 1 

The goal of the analysis is to construct the error bounds on F Y at y T . 

Assume that is normal. The cdf of X^ is written as follows noting that 
0^ * (p , g^) is a random vector. 

F i (x ; v ±t ^ (2.2) 


The best estimate of the cdf of X^ is, 

F*(x ; v ± , o i ) 



(2.3) 


Because Y is a linear function of normal X , the estimated distribution 
Y in the neighborhood of y' will be normal using the parameter estimates, 



(2.4) 
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where 


+ E a i v ± 

i-1 1 1 


‘2 

°Y 


K , .. 
r 2 2 

i-l 1 1 


(2.5) 


( 2 . 6 ) 


2.3 Properties of Y and X^ at the Design Point 

* A 

A basic property of the design point values X used to compute F^at y' is 

y’ - f(x*) 


+ A, •* x i 


i-i 


(2.7) 


where X^ is the design point associated with variable X^. 


The cdf at X is, 


* 

'X - Uj 


( 2 . 8 ) 


it it 

Shown in Fig. 2.1 is the cdf of X^ and the point X^ and F^. 

At F*. X. can be considered as a random variable denoted as X , , because it 
i* i oi 

is a function of * (u^, o ) . 


W V ■ h 1 <V 


# _1 [F*] + \i ± 


(2.9) 


Upon substituting Eq. 2.8, it follows that, 

X oi (u i’ °i ) ‘ °i X * * 

\ 

where, * 

* X - 


( 2 . 10 ) 


X 


( 2 . 11 ) 
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Note that 


ECXoi) - X 


( 2 . 12 ) 


because E^) - a ± , and E(y i > - 


Also note that X q ^ in Eq. 2.10 is a random variable because a J and 


are random variables. The pdf of X q1 Is shown In Fig. 2.1. 


2.4 The Distribution of the Response Variable at a Given F 
Define F* as the value of F corresponding to y'. 


F’ - « 


- U, 


(2.13) 


Define Y as, 
o 


a n + l a. X . 
° i-1 1 oi 


(2.14) 


*J. a i < 0 i x * + V 


i-1 


Note that from Eqs. 2.7 and 2.12, 


E(Y q ) - y' 


(2.15) 


Thus Y q is a random variable, denoting Y at F'. It is a function of 0, and 

it models or represents the error bound in Y at a given F'. The distribution 

of Y is shovn in Fig. 2.2. 

The standard deviation of Y is 

o 


o 

o 


K 

r 

L 

i-1 


(x*) 2 V( 0i ) + V( Ul ) 


1/2 


(2.16) 


And the coefficient of variation (COV) of Y is 

o 


C - 

o 


Oq/u 


o 


(2.17) 
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Fig. 2,2 The estimated distribution function for response variable Y 
and upper and lower error a, error bounds for Y 
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2 

If is normal, then will be normal and will have an x distribution. 

In general, it would be difficult to derive the distribution of Y , but as 

o 

indicated below, a "default” lognormal model can be assumed. As a general 
purpose distribution, the lognormal can be used as an approximating model 
in a variety of applications in which the exact form cannot be found. 

2.5 Upper and Lower Error Bounds on 

Option 1. The lognormal model for Y q . Assume that Y q has a lognormal 

distribution. Upper and lower error bounds on Y , denoted as y and y , and shown 

O U L 

in Fig. 2.2, can be derived as follows: Let a be the confidence level. Then, 
y^, for example, is related to o as. 


P 


Y 


< 




(2.18) 


And if Y q is lognormal, it follows that the lower error bound for Y is. 


y L * *o €XP(z l 6) 


(2.19) 


where 

'v / Y 1 

Y « yV/1 + C 
o o 

6 - /n(l + C 2 )’ 
o 


( 2 . 20 ) 

( 2 . 21 ) 


* standard normal variate at a 
probability level of (1 - a)/2 


Similarly, the upper error bound for Y q is, 

y n - ? exp(z 6) (2.22) 

u o u 

where is the standard normal variate at a probability level of (1 + a)/2. 
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Option 2. The normal model for Y q . Unfortunately, the lognormal model 

for Y has a serious limitation. The lognormal distribution is defined 
o 

only for Y > 0. If the response variable has values at zero, or in the 

"neighborhood" of zero, then the lognormal is not suitable. In general, 

this may not be a problem, but it is not unreasonable to imagine interest 

in some variable which has a zero mean. In any case, the problem can be 

avoided by using a normal model for Y. The penalty may be a loss of accuracy. 

The mean and standard deviation of Y are equal to y' and a (Eq.- 2.16) 

o o 

respectively. Then the upper and lower error bounds on Y q are, 

y L“ Z L°Y +y ’ (2>23) 

y u“ Z U a Y +y ’ (2>24) 


2.6 Translation of Error Bounds on Y to Error Bounds on 


Assuming that Y * f(£) is linear in the i neighborhood of y 1 and that all 

X,, are normal, it follows that Y will also be normal. Error bounds on Y 
i o 

can then be easily transferred to F given y f . The scheme for doing this 
is suggested in Fig. 2.3. 

The standard normal variate z defines the estimated distribution of Y, 


z 



(2.25) 


where, 

F y (y) - *(z) (2.26) 

At y' the error bounds are assumed to be parallel to as shown in Fig. 2.3. 
The slope of the line is 


dz 

dX 


1 


a 


Y 


(2.27) 
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And It is seen from Fig. 2.3 that 

(2.28) 

(2.29) 

(2.30) 

2.7 Concluding Remarks 

These first order error bounds were derived on the basis of distributional 
assumptions. It is hypothesized that these bounds are robust in that they 
provide a reasonable approximation to those bounds for the general case 
where y ■ f(X) is not linear, and X is not normal. This has yet to be proven. 

A numerical example is provided in the next chapter. 
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3.0 EXAMPLE: A NUMERICAL EXAMPLE OF FIRST ORDER ERROR BOUNDS 

Statement of the Problem 

Consider the response variable Y which is a function of R and T. 


Y « R - T 


There is uncertainty in the parameters of R and T. it is required to compute 
90% error bounds for F y , the distribution function of Y, at y' - 1 and y' - 2. 
Observations on R and T have been made . The statistics are. 


For R 
n - 20 
R - 10 


For T 
n - 20 
T - 5.0 


Thus the estimators are, 


V (U R. °R> 


S T ■ <“T, V 


"R ' 10 



Solution 

The calculations below follow the forms provided in the Summary 
The variances of the parameters are, 

V(u R ) * s R /n * (2) 2/20 - 0.20 
V(a R ) - s R /2(n - 1) - (2) 2 /38 - 0.105 
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V(u T ) - s*/n - 1/20 - 0,05 


V(a T ) - s^/2(n - 1) - 1/38- 0.026 


The design points at y' * 0 and y' ■ 1 are defined in Fig. 3.1. Thus, at 


Y - y\ 


R* - y 


R 6.8 - 10 


■ - 1.6 


^ T * “ y T 5.8-5 . „ D 

t* - 7 ■ Y — “ + 0.8 

°T 


The random variable Y given F' is denoted as Y q . The mean of Y^ is 


y o - y 


and the standard deviation of Y is, 

o 


c o - [ (r*) 2 V(o R ) + V(y R ) + (t*) 2 V(o T ) + V( Vj )] 1/2 
- [ (1 .6) 2 (1.05) + (.2) + ( .18) 2 (.026) + .05] 1/2 


o - 0.73 
o 


OPTION 1 (Lognormal model for Y ) 

o 

The C0V of Y is, 
o 

c .fo.fo . 0^3 . 0>73 
o y o y’ 1 


6 * An(l + C ) 
o 


^n(l + .73 ) * 0.654 



Y - fZ“T 12 '*• KJ ( lOj T-) 

T 

tZ-eVDCeiP VAXt\AT^l£<, t~T=- 

PC*-£- Tie; 

Luw^i f 'Sfe.4c U*. "Z n. = a- - £5 -nj> 

fZ&LuceA. Co t>r*L> 


Pe$l6U Poi/siT^ 


z.U> 


o 1 <zl> 


+-I l <7 4 


+ 7 - * 0 > 


- i.Z 


S’ 


* 7-6 


P^ 3,1 V£$i&s* Vc>wt<> CpiZfo€$1?t>rtV>lb)6i T 


m 
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The median of Y fl 



The estimated mean and standard deviation of Y, 



For 90% error bounds, (1 - o)/2 - 0.05, and (1 + a)/2 - 0.95. 

z - - 1.64 Z U “ + 1,64 

The upper and lower error bounds for Y q , 

y * y exp(z 6) * 0.80 exp (-1.64 x .654) 
L o L 

- 0.273 

* Y exp(z^6) * 0.80 exp (1.64 x .654) 

- 2.33 



Then, F' - $(-1,787) = 0.037 
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And finally the error bounds on F at y' ■ 1, 



- •(- 1 - 787 - -H^ 1 ) 

- 0.0086 


F 


U 



<- 


1.787 + 


1 - 0.273 
2.24 


■ 0.072 


These bounds are plotted in Fig. 3.2. 

OPTION 2 (Normal model for Y ) 

o 

When the normal model (Option 2) for is employed, the upper and lower 
bounds are 


- (- 1,64) (0.73) + 1. - - 0.20 

y u ■ z u °0 + y ’ 

■ (1.64) (0.73) + 1. - 2.20 
Employing the forms as above for F and F , 

L* U 

F L - 0.10 Fy - 0.105 


The Option 2 bounds also plotted in Fig. 3.2 do not agree well with the 
Option 1 bounds. Brief commentary on these differences is provided below. 


£ 



O I 7- 


1rt& Bztzptz x^otZ' Fv • 



Now compute the error bounds at y' • 2. Only those calculations which 
differ from above will be shown in the following. 

The design point (See Fig. 3.1) 


OPTION 1 


OPTION 1 


r* - - 1.2 


t* - 0.6 


o * 0.649 (Eq. ) 
o 


z • 


2-5 

2.24 


- - 1.34 


F' - 0.090 


C - 0.324 
o 

6 - 0.316 


Y - 1.90 
o 

y L - 1.11 

>’U " 3 ’ 23 
F l - 0.029 

Fy - 0.172 


y L - (- 1.64().649) + 2. - 0.935 


yy - (1.64) ( .649) + 2. - 3.06 
F l - 0.035 


0.193 


The error bounds are plotted on Fig. 3.2 
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Some Comments 

The differences in the first order error bounds reflect concerns by this 

author regarding the general use of the lognormal model for Y q . As y' -*■ 0, 

it is noted that C • and the model "blows up." 
o 

However, the lognormal may be a more accurate model for response vari- 
ables which are guaranteed to have positive values. On the other hand, 

the normal Y model avoids any mathematical difficulties in the neighbor- 
o 

hood of y' ■ 0. And because (a) the normal and lognormal error bounds 
are "reasonably close" in the region where y' > 0 and (b) the normal model 
is easier to use, it is suggested that the normal be used as a first order 
approximation to the error bounds. 
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1.0 INTRODUCTION 

1.1 Introductory Remarks 

Monte Carlo traditionally has been considered to be a "last resort" 
method for solving a probability or statistics problem because of high 
cost relative to accuracy of the results. However, in recent times a 
combination of the development of new efficient numerical techniques 
and new digital computing hardware have made Monte Carlo more attractive. 

Presented in this report are descriptions of the following Monte 
Carlo programs dedicated to probabilistic structural analysis. 

1. "Conventional" Monte Carlo 

2. Variance reduction using antithetic variates 

3. To be added later 

4. To be added later 

Provided in the following sections afe descriptions of how each method 
works as well as a comprehensive study of the performance of each. 

1 . 2 The Basic Problems 

Consider the random variable 2 as a function of the random vector 


X 

-v 


(X 


1 * 2 * 


V 


Z - h(£) 


( 1 . 1 ) 


The distribution of each is known. It is assumed that all are 
mutually independent. 

One problem of probabilistic mechanics and design is to compute a 


point probability, 


p - P[hQp - h Q ] 


( 1 . 2 ) 



For example, p could represent the probability of exceedance of a deflec- 
tion or perhaps the probability of failure. 

The second problem is the extension of the first to the construction 
of a cumulative distribution function. 

F z (z) - P[hqp - z] (1.3) 

Clearly the two problems are identical, but optimal strategies for analysis 
may differ. For example, to construct the CDF, one option would be to 
obtain point estimates of F z at selected values of z, then fit a curve 
through the points. A second option would be to construct an empirical 
distribution function from a large sample of (See Sec. 2.U). 

1 . 3 Random Samples 

The basis for Monte Carlo simulation is a standard uniform distribu- 
tion random number generator. Methods of generating uniform variates are 
generally based on recursive calculations of residues of modulus m from a 
linear transformation [ 1 ]. Most large computers have such a generator 
as a library function. 

A variety of methods can be employed to generate variates from the 
distributions. Presented in Appendix A are algorithms used for the program 
presented herein. 
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2.0 "CONVENTIONAL" MONTE CARLO 

2.1 Point Probability Estimates by Conventional Monte Carlo Using the 
Bernoulli Parameter 

Consider a function, h(£), where is a vector of random variables, 
all having known distributions. It is required to compute, 

p - P[h(j|p - h Q ] (2.1)- 

The problem can be reformulated as 

p - P[g(£) - 0] (2.2) 


where g(X), called the "performance function," is 


g(V) - h($) - h 


(2.3) 


In a direct Monte Carlo scheme, a sequence of K random vectors, 

X^, can be sampled, and in turn, a sequence of g^; i ■ 1, K computed. Define 


1 if g. - 0 

0 if gi > 0 


(2.4) 


Thus, Y^ has a Bernoulli distribution 

P(Y i - 1) ■ 


(2.5) 


P(Y i - 0) - 1 - p 


where the Bernoulli parameter p is the same p as in Eq. 2.1. 
The maximum likelihood estimate (MLE) of p is [ 5 ] , 

- 1 K 
i«l 


P 


1 

K 


( 2 . 6 ) 



But Ty. is just the total number of g, - 0, denoted as N . 
L i i o 

just the fraction of the g^'s less than zero 


P 



Thus, p is 


(2.7) 


A flow diagram of conventional Monte Carlo is given in Fig. 2.1. 

A listing of a computer program for conventional Monte Carlo employing 
the Bernoulli parameter is provided in Appendix B and an example of the 
output is shown in Fig. 2.2. 


2.2 Confidence Intervals on the Bernoulli Parameter, p 

The MLE of p is p. Because of sampling error, p is only an estimate, 
and the key question is how close is p to p. Confidence intervals are 
described below. Note that these confidence intervals refer to 
sampling error of the Monte Carlo process, not uncertainties associated 
with the parameters of X^. 

Consider p, 

1 K 

P - 7 [ Y (2.8) 

K i«l 1 

The mean and variance of p are [ 5] 

E(p) - p (2.9) 

V(P) * £ ~ ( ' 1 ~ k ~ ~ ~ E) ( 2 . 10 ) 

By the central limit theorem, p will approach a normal distribution as 
K -► ®. Confidence intervals for p are constructed using normal distribution 


mathematics, 




Fig. 2.1 Flow diagram of conventional Monte Carlo 












MONTH CARLO SOLUTION 


LIMIT STATE FUNCTION s R=S 


SAMPLE SIZE, K« 100 
NUMBER OF RANDOM VARIABLES, N- 2 

RANDOM VARIABLES 

VARIABLE DISTRIBUTION MEAN 

R WEI BULL . 20040E+D2 

S EVD . 1 0000B+02 

Note that Y is the same as g(X); 
these are the statistics on the 
limit state function. 


This is p 

T 

NUMBER OF NEE Y VALUES' 0. PERCENT OF TRIALS' .000000 

Fig. 2,2 Oucput of conventional Monte Carlo program. (No sorting requested) 
Performance function; g(R,S) * R - S 


STATISTICS OF Y i 


MEAN 
STD DEV 
MEDIAN 
COV 


. 1 00 1 BE +02 
. 27499E+0 1 
. 96606E+0 1 
. 27450E+0O 


STD DEV 
. 20000E+01 
. 20000E+O1 
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/ p(l - p) < < '' , 

P - \!2 / - P " P + Z a/2 


'(1 ~ P) 
K 


( 2 . 11 ) 


where p is substituted for p in the variance. The probability that p will 

be bounded by the lower and uppper limit is 1-a, where a is the confidence 

coefficient, z is the standard normal variate corresponding to a/2. 
a/2 

Commonly used values 


a 

Z a/2 

.10 

1.64 

.05 

1.96 

.01 

2.58 


The confidence interval of Eq.2.11 relies on the central limit theorem 
and must be considered as only an approximation for finite K. In general, 
the approximation is considered "valid" if Kp > 5 [ 5]. 

Eq. 2.11 can be written as, 1 

p (1 - Y) - p - p(l + Y) (2.12) 


where, 



(2.13) 


Eq. 2.13 is displayed in Figs. 2.3 and 2.4 for 90% and 95% confidence 
intervals respectively. These figures show the sample size requirements 
for confidence intervals of a given width and level. For example, if the 

_3 

point probability is expected to be about 10 , and it is required to have 

p within - 10% of p with a confidence of 90%, then it is necessary to have 
a sample of size K > 200,000. 






ction of sample size and p. 
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2.3 Computer CPU Time on the CYBER 175 

The conventional Monte Carlo program of Appendix B was exercised on 
several problems using all five of the available distributions. CPU time 
was recorded for each program. It is assumed that this conventional Monte 
Carlo program will provide an upper bound to CPU time relative to other, 
and more efficient, Monte Carlo schemes. The CYBER 175 is the mainframe 
computer at the University of Arizona, and all results relate to this machine. 

Recorded CPU time for several examples was consistent. Compilation and 
loading time for all cases are shown in Table 2.1. These are average values, 
but there was little variation. 

Execution CPU time essentially depends only upon the number of variables 
and not on distributional forms or performance functions. Fig. 2.5 illustrates 
the CPU execution time per variate as a function of sample size K. Total CPU 
time is obtained by adding compilation and loading time to execution time. 

A sample program was run on both the CYBER 175 and the VAX 11/780 for 
a time comparison. The results shown in Table 2.2, reaffirm the fact that 
the VAX is too slow for production Monte Carlo. 

To get an idea of computer charges for running Monte Carlo, Fig. 2.6 
is provided. This is the commercial rate of the UA CYBER 175 for low priority 


jobs . 
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Table 2.1 

Compilation and Loading CPU Time for Conventional 
Monte Carlo on CYBER 175 Program 



CPU Time (sec) 

Compile 

1.0 

Load 

0.25 


Table 2.2 

Comparison of CPU time Between CYBER 175 .and VAX 11/780 

* 

for one Example Problem 

Time (sec) 



CYBER 175 

VAX 11/780 

Compile 

1.0 

14 

Link 

0.25 

5 

Execution 

7.5 

30 




TOTAL 

8.75 

49.0 


There were 2 variables; K * 30,000. 
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COST , D 



Fig. 2.6 Cose in dollars ($), D, as a function of time for 
the UA CYBER 175; lowest priority. 
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2.4 Comparison of Monte Carlo to Wu/FPI 

Computational efficiency was the motivation for the development of the 
Wu/FPI program. It is generally known that Monte Carlo is Inefficient 
relative to a fast probability integration method. An attempt is made here 
to quantify differences in computer time between conventional Monte Carlo 
and Wu/FPI. Because the cost of conventional Monte Carlo depends upon the 
accuracy and probability level required, a general direct comparison can't 
be made. However, an example presented in the following clearly demonstrates 
the high cost of Monte Carlo. 

Suppose that it is required to provide a Monte Carlo solution such 
that the 90% Cl for p is within * 10% of p. The CPU execution time for the 
CYBER 175 can be computed from Figs. 2.3 and 2.5 for a given probability 
level, p. This CPU time is shown in Fig. 2.7 as a function of the number 
of variables in g(£) for p - 10 and 10 . At these levels Monte Carlo is 
two to three orders of magnitude more expensive than FPI. And the FPI 
solution is likely to be more accurate. Moreover, for smaller tail proba- 
bilities FPI gets no more expensive while Monte Carlo will break the bank. 

? . 5 Estimating the CDF of a Random Function 
2.5.1 The Empirical CDF 

Conventional Monte Carlo provides capability for estimating the complete 
distribution function of a function of random variables . Define the random 
variable 2, as a function of the random vector 


Z - Z(£) 


(2.14) 


E- 15 



15 
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A random sample of i ■ 1, K is used Co generate a random sample of 
1 ■ 1, K. In curn, an empirical distribution function of Z can be 
constructed using methods of probability plotting. The empirical CDF, 


denoted as F^, will be an estimate of the CDF of Z, F^Cz). 

Various forms of F^ have been proposed [3, 4, 6 ] . The values of 

F^ below correspond to Z^ where Z^^ is the ith smallest value of the 


random 


vector £. Thus, F^ = F^(Z.^j). 


1. Hazen; F^ * 


i - 1/2 


2. Gumbel; F i ■ 

3. Median ranks, F. » - ^7 

1 n + u.h 


Through prior experience on extensive Monte Carlo simulation, this author 
has found that the Hazen formula consistently provides "good estimates" 


2.5.2 The Sort Routine 

To construct the empirical CDF it is required to sort the random 
sample £ to obtain an ordered sample £ . Let Z... denote the ith smallest 


value. 


The routine used in this Monte Carlo code is program QUICKSORT which 
is considered to be the fastest available [7]. A description of QUICKSORT is 
given in Appendix C . The Fortran statements for this code are provided 
in the program listing in Appendix B . 

CPU time requirements for the sort routine can be relatively large for 
large samples. Fig. 2.8 shown CPU execution times as a function of the 
size of the £ vector. 




2.5.3 An Example. 

Shown in Fig. 2.9 is a cable of the sorted vector and the corres- 

ponding for the example of Fig. 2.1. This is the data required for 
plotting. The empirical CDF of Fig. 2. 10 was done by hand, but in general 
such graphs can be automated using a computer graphics package. 
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SORTED VALUES OF Z AND THE EMPIRICAL CDF 



m 

l 

. 32159E+01 


m 

6 

. 48457E+01 


m 

1 1 

. 59944E+01 


m 

16 

. 69827E+01 


m 

21 

. 76156E+01 


s 

26 

. 87304E+01 


m 

31 

. 90619E+01 


= 

36 

. 92816E+01 


m 

41 

. 95862E+01 


s 

46 

. 10054E+02 


s 

51 

. 10376E+02 


m 

56 

. 10712E+02 


m 

61 

. 10856E+02 


m 

66 

. 1 1 191E+02 


K 

71 

. 1 17308+02 


S 

76 

. 12122E+02 


m 

81 

. 126678+02 


m 

86 

. 1 2893E+02 


m 

91 

. 132738+02 


m 

96 

. 13943E+02 



. S 

l 

. 500008-02 



6 

.550008-01 


■ 

11 

. 105008+00 


* 

16 

. 155008+00 


s 

21 

. 205008+00 


* 

26 

.255008+00 



31 

. 305008+00 


= 

36 

.355008+00 



41 

.405008+00 


* 

46 

. 45500E+00 


= 

51 

. 50500E+00 


= 

56 

.555008+00 



61 

. 60500E+00 


* 

66 

. 65500E+00 


* 

71 

. 70500E+00 


* 

76 

. 75500E+00 


St 

81 

. 80500E+00 


* 

86 

. 8S500E+00 


* 

91 

.905008+00 


« 

96 

. 95500E+00 


. 40876E+01 • 4283 IE ■*■01 
. 48984E+01 . 50586E+01 
. 60426E+0 1 . 66202E+0 1 
. 70597E+01 . 70685E+0 1 
. 79653E+01 . 8386 IE +01 
. 87709E+0 1 . B7964E+0 1 
. 9097 1E+01 . 91454E+01 
. 92823E+01 . 93259E+B1 
. 95993E+0 1 . 96380E+0 1 
. 10115E+02 . 10137E+02 
. 1 058 1 E+02 .1 0607E+02 
. 10771E+02 . 10773E+02 
. 10874E+02 . 10958E+02 
. 1 1 2468+02 . 1 1 344E+02 
. 1 1 760E+02 . 1 1 B02E+02 
. 12140E+02 . 122B4E+02 
. 12803E+02 . 12844E+02 
. 12963E+02 . 13042E+02 
.132978+02 . 13361 E+02 
. 14797E+02 . 14983E+02 


. 1 5000E-0 1 . 25000E-0 1 
. 65000E-01 . 75U00E-01 
. 1 15008+00 . 12500E+00 
. 165008+00 . 17500E+00 
. 21500E+00 . 22L00E+00 
. 26500E+00 . 27500E+00 
. 31500E+00 . 32500E+00 
. 36500E+00 . 375008+00 
. 41500E+00 . 42500E+00 
. 465008+00 . 47500E+00 
. 5 1 500E+00 . 52500E+00 
. 56500E+00 . 57500E+00 
. 61500E+00 . 62500E+00 
. 66500E+00 .675008+00 
. 71500E+00 . 72500E+00 
. 765 00E+00 . 77500E+00 
.81500E+00 .825008+00 
. 865008+00 . B7500E+00 
. 915008+00 . 92500E+00 
. 96S00E+00 . 97500E+00 


. 44764E+01 .45626E+01 
.5615BE+01 . 59102E+01 
. 68500E+01 .69210E+01 
.70780E+01 . 710048+01 
.84534E+01 . 84720E+01 
. 888508 + 01 . 89137E+01 
. 92372E+01 . 925578+01 
. 95770E+0 1 . 95829E+01 
. 9S157E+01 . 98782E+01 
. 10256E+02 . 10370E+02 
. 1B631E+02 . 10644E+02 
.1079 1 E+02 . 1 0846E+02 
. 1U25E+02 . 1U62E+02 
. 1 14098+02 . 1 16168+02 
. 11912E+02 . 1 1933E+02 
.1241 3E+02 . 1 2573E+02 
. 1 2B67E+02 . 128738*02 
. 13131E+02 . 13142E+02 
. 1363BE+02 . 13709E+02 
. 15123E+02 . 15305E+02 


.350008-01 .450008-01 
. 85000E-01 .950008-01 
. 1 3500E+00 . 1 450UE+00 
. 18500E+00 . 195008 -*-00 
. 23500E+00 .245008-1-00 
. 285008+00 . 29500E+00 
. 335008+00 . 345008+00 
. 38500E+00 .395008+00 
. 43500E+00 . 44500E+00 
.485008+00 .495008+00 
.535008+00 . 54500E+00 
•585008+00 .595008+00 
.635008+00 . 64500E+00 
. 685008+00 . 695008+00 
.735008+00 .745008+00 
. 78500E+00 . 795008-*00 
.835008+00 . 84500E+00 
. 88500E+00 . 69500E+00 
.935008+00 .945008+00 
.985008+00 .99500E+00 


Fig. 2.9 Sorted and corresponding empirical CDF for the example of Fig. 2.1 
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3.0 THE VARIANCE REDUCTION METHOD 

3.1 Preliminary Remarks 

The variance of Monte Carlo estimators can be reduced, relative to 
straightforward sampling of Chapt. 2.0, by appropriate operations with 
negatively correlated samples. Ang and Tang [ 1] present several examples 
which demonstrate dramatic improvements in efficiency realized by variance • 
reduction methods. 

A variance reduction computer program, tailored for structural 
mechanics analysis by providing point probability estimates of functions of 
random variables has been developed. The listing is given in Appendix D. 

To assess performance, the program has been exercised on several examples. 
Results presented in Section 3.6 show dramatic improvement of variance 

reduction over conventional Monte Carlo in some cases. In other cases, 

t 

the improvement is only modest . Some general conclusions are presented 
in Section 3.7. For the most part however, for a given problem it is dif- 
ficult to predict how much improvement one can expect with variance reduc- 
t ion . 


3.2 The Essence of Variance Reduction 
The goal of analysis is to estimate 

p - P[hqp < h o ) (3.1) 

Suppose p and p' are two unbiased estimates of p. (The method for obtaining 
a point estimate of p is described in Sec. 3.4 below.) The two estimators 
may be combined to form another estimator 

P * ^<P + P f ) 


(3.2) 
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The expected value of p £ is, 

E(p) - j[E(p) + E(p') ] - p 

which means that p is an unbiased estimator. 

The corresponding variance is 

V(p) - ^lV(p) + V(p') + 2 Cov (p, p’)] 


(3.3) 


(3.4) 


If p and p' are statistically independent, for example, based on two separate 
and independent sets of random numbers, 

V(P) - J [V(p) + V(p')J (3.6) 


Thus, the accuracy of the estimator p can be improved over that of the 
independent case if p and p' are negatively correlated. Ang and Tang cite 
several examples (no structural analysis) where variance reduction can 
provide a dramatic improvement in efficiency of probability estimation [1J. 

An estimate of p is obtained by several samples, p^; i » 1,K. 


P 


E 


i K 

i? 

K i-i 1 


(3.7) 


all p^ are independent. Note that p^ will approach normality as K + • 
as a consequence of the central limit theorem. 

The mean and variance of p E are. 


E(P E ) - P 

V(p_) - o 2 /K 
E p 

2 

where o is estimated as, 

P 


2 

s 

P 


1 

K-l 


i-1 


(Pi > 


P E } 


2 


(3.8) 

(3.9) 


(3.10) 
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3.3 How to Obtain Negatively Correlated Samples 

Suppose chat Che uniformly distributed variate u^ Is used to generate 

a number from a given distribution (See Appendix A). Then the uniform 

variate u' ■ 1 - u will produce an x' such that x and x' will be negatively 
i i 1 1 i 

correlated. The u| are called "antithetic" variates. 

And in general, if u^, u^ , . . . u r is used to generate p, and 1 - u^, 

1 - u^, . . . 1 - u is used to generate p', then p and p' will be nega- 
tively correlated. 

Such a procedure works well when the integral transform is used, e.g., 
Weibull, EVD. One uniform variate is used to generate. one x^. But 
where Box-Muller is used to generate normal' variates , two u^ are chosen 
(See Appendix A). While the resulting x^ and x^ will be negatively correlated, 
the correlation coefficient will not be -1.0. An improvement can be made 
by choosing x!^ as a "mirror image" of x^ in the distributions. This can 
be done by 

x^ * 2u - x^ (3.11) 

where y is the mean of X. 

3.4 How to Obtain Point Probability Estimates 
3.4.1 The Two Variable Case 

The structural reliability problem in which p is the probability of 
failure will be used to illustrate how p and p' are obcained. Consider 
the design case where the two variables are R (strength) and S (stress). 
Estimate p, where 


p - P[R - S - 0] 


(3.12) 
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Both R and S are random variables whose density functions are shown 
in Fig. 3.1. First S, having been identified as the variable having the 
largest variance, is the "reference." A random variate is sampled from 
the other factor, R. An estimate of p is 

p £ - P(S > K ± ) 

- 1 - F s (R i ) (3.13) 

where Fg is the CDF of S. 

It should now be apparent why sampling is done on the smallest vari- 
ance term, p is a "good" estimate of p if the distribution is narrow, and 
is exact as + 0. 

Now the antithetic variate R^ is sampled as described above. Because 
it is negatively correlated to R^, its position relative to R^ will be as 
shown in Fig. 3.2. Then, 

pj - P(S > R^) 

■ 1 - F (R!) (3.14) 

u X 

and the ith estimate of p is 

h ■ I <fi + ?;> (3 - 15) 

As a second example, consider again the case where R and S are the basic 

variables, but now where o D < o_. In this case, R would be the reference 

K 5 

variable. Random points S 4 and the antithetic variate S! are sampled from 
S. The estimates now are, 

^ - W 
pi - V s i> 


(3.16) 


1 — - P A - P(s > R 1 ) 

Fig. 3.1 Estimate of p using one point sampled from the minimum 
variance variable. 


P(s > R t ) 


PCS > R’) 


Fig. 3 


2 Estimates of p using a point R^ sampled from R and the 
antithetic variate of R,, denoted as R’ 
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Thus, it is seen that the variable type (stress or strength) must be identi- 
fied to obtain the proper form for computing estimates. 

Fig. 3.2 shows why negatively correlated variables tend to provide 
good estimates. Being on both sides of a distribution, and Rj^ combine to 
produce an "average" estimate of p. 


3.4.2 The General Case 

In general, the performance function, g(Jt) ■ h(^) - h^ is a non-linear 
function of several variables. The method of obtaining a point estimate of 
p is an extension of the scheme for two variables. 

The reference variable is defined, not as the one having the maximum 
variance, but rather the one having the maximum impact. For example, if 

g » 5R - S (3.17) 

and o * o_/2, clearly the random variable, iR. * 5R will have a larger vari- 
R S ^ 

ance than S. Thus, we say that R is the maximum impact variable. 

In general, the maximum impact variable can be found by estimating 
9g/3X i for each The maximum impact variable, denoted as X^, is that 

X^ for which | 3 g / 3X ^ | is the largest. 

The sign of 3g/3X i identifies variable type; stress if (+) and strength 
if (-). As indicated above, the "type" of must be known to choose the 

appropriate form for estimating p (e.g., Eqs. 3.13 and 3.16). 

The estimates p and p f proceed as follows. Sample all variables but 

X^. Let g(J£) * 0, and solve for (this is done by the secant method 

in the program) . 


x 


M 




(3.18) 


where x is the vector of sampled Q minus X^. 
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The estimate of p is. 


Fx^(x m ) if Xj^ is a strength variable 


1 - F„ (Xj^) if ^ is a stress variable 


(3.19) 


V 


To obtain p', the antithetic vector of ^ is used in Eq. 3.19, 


3.5 Confidence Intervals on p 

Noting that p E is normally distributed, approximate 1 -a confidence inter- 
vals on p can be constructed as [ 5 ], 


p e --7=-^<p< p e + ^-^ 

/k rr 


(3.20) 


or, 


p £ (l - y) < p < p £ (l + y) 


(3.21) 


where , 


z ^/2 * standard normal variate (absolute value) at 
probability level a/2. 


Z a/2 C P 

v r K~ 


(3.22) 


C p“ S P /P E 


(3.23) 


The UA variance reduction program chooses K to produce a specific 
confidence interval. For example, if you want to sample until the 952 
confidence intervals are - 102 of p E , 


Y - 0.10 


V2 * l - 64 


0.24) 
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and solving Eq. 3.22 for K, 

C 2 

K > — E — . 269 C 2 (3.25) 

Y P 

To find C . an initial sample of K - 1000 is chosen and an estimate 
P 

of C is obtained. Then if K < 1000 in Eq. 3.16, the process is terminated 
P 

with narrower confidence intervals than requested. If K > 1000, the. program 
will continue to sample to that value. 

3.6 The Variance Reduction Monte Carlo Program 

A flow diagram which outlines the logic of the variance reduction 
program is provided in Fig. 3.3. Sample output of the program is shown 
in Fig. 3.4 with some commentary. 

Two versions of the program have been developed. An interactive version 
(IVARED) runs on the IBM PC/XT. Program VARED runs on the VAX or CYBER 175. 

A listing of VARED is given in Appendix D. 

3. 7 Examples of the Performance of VARED 

Twelve examples of the use of VARED to produce point probability estimate 
are provided in Tables 3.1 through 3*12. Point estimates by VARED are compared 
to the exact solution (closed form or P0FAIL) if available. The exact 
solution, provided by program P0FAIL, is employed for performance functions 
involving two variables. For larger problems, Wu/FPI is used. For the 
VARED solutions, 95X confidence intervals (a - 5%) are specified along 
with y “ 0.10. 

To compare variance reduction with conventional Monte Carlo, sample 
size requirements and CPU time for the latter are extracted from Figs. 2.4 and 
2.5 and are presented in the tables. 





( 1 ) 


( 2 ) 

(3) 

(4) 

(5) 

( 6 ) 
(7) 



( 8 ) 


? 


Repeat steps (3) through (7); 1 • 1,K 



Fig. 3.3 An outline of the variance reduction Monte Carlo program 
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Fig. 3.4 An example of the output of the variance reduction Monte Carlo 
Program with commentary 

MONTE CARLO SOLUTION 


LIMIT STATE FUNCTION s G*R-DSQRT (300. *P**2+1 . 92*T*#2> 


SAMPLE SIZE - 1000 

NUMBER OF RANDOM VARIABLES « 3 

CONFIDENCE INTERVAL » 95.00 */. 

GAMMA - .10 

MAX. IMPACT VARIABLE ■ X< 1) 
VARIABLE TYPE IS STRENGTH 



1 , 


VARIABLE 

R 

P 

T 


RANDOM VARIABLES 
DISTRIBUTION MEAN 

WE I BULL . 48000E+02 

LOG . 9B700E+00 

EVD . 20000E+02 


STD DEV 
. 30000E+0 1 
. 16000E+00 
. 20000E+01 


ESTIMATE OF P « . 16043E-02 

95.00 V. CONFIDENCE INTERVALS ARE 
PL « . 1 1 725E-02 PU « 


This is the first estimate of p 


20360E-02 


Note that 95% confidence 
intervals exceed - 10%. 
Thus, a larger K is 
required. (See below) 


STATISTICS OF P i 


MEAN 


1 6043E-02 
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STD DEV ■ 
MEDIAN « 
COV - 


. 69662E-02 
. 36004E-03 
.43422E+01 


K FOR BAMMA « . 10 IS 7244 

ESTIMATE OF P ■ . 1B030E-02 

95.00 7. CONFIDENCE INTERVALS ARE 


Based on the first sample of K * IOC 
this is the total K required for th< 
desired confidence intervals. K is 
computed from Eq. 3.25 which requirf 
0^. This is why the first sample oj 
1000 is taken. 


PL « . 15509E-02 PU - .20550E-02 

Note that the confidence intervals do not quit 

meet the specifications. This is because the 

original estimate of C “4.34 was small relat 

P 

to the improved estimate of C * 6.24 

P 

STATISTICS OF P * 

MEAN « ' . 1B34BE-02 

STD DEV = . 1 1456E-01 

MEDIAN » . 290 1 7E-03 

COV * . 62436E+01 


DO YOU HAVE ANOTHER DATA SET ?<Y/N> 

Note: The size of the sample required K depends upon (Eq. 3.25). 

In this problem C is relatively large implying that a relatively 
P 

large K is required. This same problem is presented in Table 3.7. 
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Table 3.1 Example of the Performance of a Variance Reduction Monte Carlo 
Program; EXAMPLE 1 

DEMONSTRATING THE PERFORMANCE OF THE UA VARIANCE REDUCTION MONTE CARLO PROGRAM 

EXAMPLE 1 

PERFORMANCE FUNCTION: ^ = R”S 


Variable 

Type 

Mean/Median 

Std. Dev./ COV 

R 

N. 

5 0- 

-5. 

c 

^ J 

N 

2 0. 

I 2. 
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RESULTS: 



Probability 
of Failure 

Total 

'fcPUrTime( b > 

Sample 
Size, K' C/ 


1 ■ 05 1 E '2 

x 

x 

Monte Carlo 

Variance 

Reduction(^) 

i i 18 E-Z 

2 04 - 

1 &0 



n.z 

5 


Notes : 

(a) Exact value using POFAIL if two variables. If more than two, 
VJu/FPI is used; the exact should be within 5Z of this value. 

(b) CYBER 175 

(c) The number of p t for variance reduction and the number of Z i for 
conventional. The values are not directly comparable. 

(d) 952 confidence intervals within - 10Z of P E 

(e) Same confidence interval as variance reduction. 








Table 3.2 Example of the Performance of a Variance Reduction Monte Carlo 
Program; EXAMPLE 2 

DEMONSTRATING THE PERFORMANCE OF THE UA VARIANCE REDUCTION MONTE CARLO PROGRAM 
EXAMPLE 2m 

PERFORMANCE FUNCTION; R~S 


Variable 

Type 

Mean/Median 

Std. Dev./ COV 

R 

LN 

_ *50. * 

■' 0-2 * 

S 

LN 

20. x 

0Z -* 
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RESULTS: 



Probability 
of Failure 

Total | 

CPU- Time < b > 

Sample 
Size, K (c) 

Exact ^ 
Wu/FPI 

53¥\ E -if 


x 

Monte Carlo 

Variance 

Reduction^) 

5-072E-4- 

13-76 

1 1 589 

Monte Carlo 
Conventional 
(Bernoulli 
parameter) 

>< 

. 238-9 

1-122 E.6 


Notes : 

(a) Exact value using POFAIL if two variables. If more than two, 
Wu/FPI is used; the exact should be within 52 of this value. 

(b) CYBER 175 

(c) The number of for variance reduction and the number of 2^ for 
conventional. The values are not directly comparable. 

(d) 952 confidence intervals within - 102 of p^ 

(e) Same confidence interval as variance reduction. 
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Table 3.3 Example of the Performance of a Variance Reduction Monte Carlo 
Program; EXAMPLE 3 

DEMONSTRATING THE PERFORMANCE OF THE UA VARIANCE REDUCTION MONTE CARLO PROGRAM 

EXAMPLE 

PERFORMANCE FUNCTION: ^ = R 


Variable 

Type 

* 

Mean/Median 

Std. Dev./ COV 

ft 

WEI. 

4.5 

0. 45 

S 

FRE. 

3-0 

0.30 















, 



RESULTS: 


V 



Sample 
Size, K (c) 

Exact 

Wu/FPI 

i.o933£-2 

x 

x 

Monte Carlo 

Variance 

Reduccion(d) 



>7 3 7 

Monte Carlo 
Conventional 
(Bernoulli 
parameter) 





Notes : 

(a) Exact value using POFAIL if two variables. If more than two, 
Wu/FPI is used; the exact should be within 5% of this value. 

(b) CYBER 175 

(c) The number of p^ for variance reduction and the number of for 
conventional. The values are not directly comparable. 

(d) 95% confidence intervals within - 10% of p £ 

(e) Same confidence interval as variance reduction. 






Table 3.4 Example of the Performance of a Variance Reduction Monte Carlo 
Program; EXAMPLE 4 

DEMONSTRATING THE PERFORMANCE OF THE UA VARIANCE REDUCTION MONTE CARLO PROGRAM 
EXAMPLE 4- 

PERFORMANCE FUNCTION: ^ js 


Variable 

Type 

Mean/Median* 

Std. Dev./ COV* 

R 

VJe i . 

>9- 


S 

P Y«. 


0- b 


















RESULTS: 



Probability 
of Failure 

Total 

CPU-Time( b > 


Exact ^ 
Wu/FPI 

<f.V)3 £-2 

>c 

HH 

Monte Carlo 

Variance 

Reduction^) 

U.oSil £-2 

}^4 $ 


Monte Carlo 
Conventional 
(Bernoulli 
parameter) 

f>< 

M8<) 

1 


Notes : 

(a) Exact value using POFAIL if two variables. If more than two, 
Wu/FPI is used; the exact should be within 5% of this value. 

(b> CYBER 175 

(c) The number of for variance reduction and the number of for 
conventional. The values are not directly comparable. 

(d) 95 7, confidence intervals within - 10Z of p^ 

(e) Same confidence interval as variance reduction. 


















Table 3.5 Example of che Performance of a Variance Reduction 1*^ c«lo 
Program; EXAMPLE 5 

demonstrating the performance of the ua variance reduction monte carlo program 

EXAMPLE 


PERFORMANCE FUNCTION: J r R'S 


Variable 

Type 

Mean/Median* 

Scd. Dev./ COV* 

R 

■A! e,i 

20 . 

2.0 

S 

EVP 

1 o. 

2.0 














— 




RESULTS: 



Probability 
of Failure 

Total 

CPU-Tipe^) 

Sample 
Size, K< c ) 

Exact ^ 
Wu/FPI 

- pr^r-^fT.i 



Monte Carlo 

Variance 

Reduction^) 

G (-73 E-? 

ic?-86l 

1 II ?fc2 

Monte Carlo 
Conventional 
(Bernoulli 
parameter) ( e ) 

X 

36- 151 

1 O 


Notes : 

(a) Exact value using POFAIL if two variables. If more chan two, 
Wu/FPI is used; the exact should be within 52 of this value. 

(b) CYBER 175 

(c) The number of for variance reduction and the number of for 
conventional. The values are not directly comparable. 

(d) 95 2 confidence intervals within - 102 of p 

r E 

(e) Same confidence interval 


as variance reduction. 
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Table 3.6 


Example of the Performance of a Variance Reduction Monte Carlo 
Program; EXAMPLE 6 


DOiONSTRATING THE PERFORMANCE OF THE UA VARIANCE 
EXAMPLE fi 


REDUCTION MONTE CARLO PROGRAM 


PERFORMANCE 

Variable 

FUNCTION: ^ r 

Type 

BTsl ,s - 

Mean/Median 

Std. Dev./ COV* 


LN 

1 * 0* 

0- 3* 



WEI 

43 E9 

0. 

_ i 

LN 

0-9* 

— v Y 

0.2 5* 














RESULTS: 



Probability 
of Failure 

Total 

CPU-Time^) 

Sample 
Size, K< c > 

Exact 00 

Wu/FPI 

1-901 E-3 

7>~<f 


no nee ua r l o 

Variance 

Reduction(d) 

I-19E8E-3 


1 14-37 

Monte Carlo 
Conventional 
(Bernoulli 
parameter) ( e ) 


b8-3Gi6 

13952 6 


Notes : 

(a) Exacc value using POFAIL if two variables. If more than two, 
Wu/FPI is used; the exact should be within 5Z of this value. 

(b) CYBER 175 

(c) The number of for variance reduction and the number of z for 
conventional. The values are not directly comparable. 

(d) 95 7. confidence intervals within - 10Z of d 

r E 

(e) Same confidence interval as variance reduction. 
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Table 3.7 Example of the Performance of a Variance Reduction Monte Carlo 
Program; EXAMPLE 7 

DEMONSTRATING THE PERFORMANCE OF THE UA VARIANCE REDUCTION MONTE CARLO PROGRAM 

EXAMPLE *"7. 


PERFORMANCE ! 
Variable 

FUNCTION: ^ - 
Type 

p a i \.<\2 ■ 
Mean/Median 

P 1 

Std. Dev./ COV* 

R 

Wei 

48-0 

3-0 

P 

LNJ 


0-1 6* 

T 


20-0 

2.0 














RESULTS: 



Probability 
of Failure 

Total 

CPU-TimeC b ) 

Sample 
Size, K (c > 

Exact 

Wu/FPI 

0-00 16 



Monte Carlo 

Variance 

Reduction^) 

0 0018208 

16-575 

12734- 

Monte Carlo 
Conventional 
(Bernoulli 
parameter) ( e ) 


74-4-166 

211349 


Notes: 

(a) Exact value using POFAIL if two variables. If more than two, 
Wu/FPI is used; the exact should be within 5Z of this value. 

(b) CYBER 175 

(c) The number of for variance reduction and the number of 2^ for 
conventional. The values are not directly comparable. 

(d) 95 7. confidence intervals within - 10Z of 

(e) Same confidence interval as variance reduction. 


Table 3.8 Example of the Performance of a Variance Reduction Monte Carlo 
Program; EXAMPLE 8 


DEMONSTRATING THE PERFORMANCE OF 


EXAMPLE 



THE UA VARIANCE REDUCTION MONTE CARLO PROGRAM 


PERFORMANCE 

FUNCTION : *3 — 

A- 10000 . , £ ff 

i I - 4 f? 

Variable 

Type 

1 b(Y*< 

Mean/Median* 

ro)"* 11 

Scd. Dev./ COV 

A 

IN. 



f rr 

N. 

0-7 

o • o 7 

G 

LN. 

0 - 227 * 

0 ' o-* 

Y 

r ln 

LO* 

1 0 , 15 * 

AG 

EVP. 

o • ooo5 

v-ocooft 

H 

LN 

1 

0 ^* 


RESULTS: 



Probability 
of Failure 

Total 

CPU-Time( b ) 

Sample 

Size, 

Exact (a) 
Wu/FPI 

1-002 B-2 


x 

Monte Carlo 

Variance 

Reduction(^) 

1-881 4- E-j 

i4-8^ 

M-ol 

Monte Carlo 
Conventional 
(Bernoulli 
parameter) ( € ) 



19810 


Notes : 

(a) Exact value using POFAIL if two variables. If more than two, 
Wu/FPI Is used; the exact should be within 5Z of this value. 

(b) CYBER 175 

(c) The number of p^ for variance reduction and the number of Z for 
conventional. The values are not directly comparable. 

(d) 95 Z confidence intervals within - 10X of d 

r E 

(e) Same confidence interval as variance reduction. 






















Table 3.9 Example of the Performance of a Variance Reduction Monte Carlo 
Program; EXAMPLE 9a 

DEMONSTRATING THE PERFORMANCE OF THE UA VARIANCE REDUCTION MONTE CARLO PROGRAM 

EXAMPLE _ 3 a 

PERFORMANCE FUNCTION: z 


Variable 

Type 

Mean/Median* 

Std. Dev./ COV* 

R 

LN 

20 0* 

■ 0^* 

S 

LN 

1 0.0* 





* 














RESULTS: 



Probability 
of Failure 

Total 

CPU-Time( b ) 

Sample 
Size, K< c > 

Exact ^ 
Wu/FPI 

(?• hb'^2 


1 

Monte Carlo 

Variance 

Reduction^) 

k IH596-? 

4~'l5 

30? 1 

Monte Carlo 
Conventional 

jf 

i 5-172 4- 


(Bernoulli 
parameter) ( e ) 

Bfl 


Notes : 

(a) Exact value using POFAIL if two variables. If more than two, 
Wu/FPI is used; the exact should be within 5X of this value. 

(b) CYBER 175 

(c) The number of p^ for variance reduction and the number of for 
conventional. The values are not directly comparable. 

(d) 95Z confidence intervals within - 10X of p 

w 

(e) Same confidence interval as variance reduction. 









Table 3.10 Example of the Performance of a Variance Reduction Monte Carlo 
Program; EXAMPLE ?b 

DEMONSTRATING THE PERFORMANCE OF THE UA VARIANCE REDUCTION MONTE CARLO PROGRAM 
EXAMPLE [2 

PERFORMANCE FUNCTION: ^ = R"S 


Variable 

Type 

Mean/Median* 

Std. Dev./ COV* 

R 

LN 

22 - 5 y 

Ov2 * 

5 

LN 

1 0 . 0 * 

0 , z y 















- 



RESULTS: 



Probability 
of Failure 

Total 
CPU- Time 

Sample 
Sire, K< c > 

Exact' fa> 
Wu/FPI 

1-89336 E-3 



Monte Carlo 

Variance 

Reduction^) 

l-'WHE-? 

9- o') 5 

bob? 


x 

Sl-W 

>| 877 b 


Notes: 

(a) Exact value using POFAIL if two variables. If more chan two, 
Vu/FPI is used; the exact should be vichin 51 of this value. 

(b) CYBER 175 

(c) The number of p^ for variance reduction and the number of for 
conventional. The values are not directly comparable. 

(d) 957 . confidence intervals within - 10X of p £ 

(e) Same confidence interval as variance reduction. 
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Table 3.11 Example of the Performance of a Variance Reduction Monte Carlo 
Program; EXAMPLE 9c 

DEMONSTRATING THE PERFORMANCE OF THE UA VARIANCE REDUCTION MONTE CARLO PROGRAM 
EXAMPLE q C 

PERFORMANCE FUNCTION : ^ r R. 5 


Variable 

Type 

Mean/Median* 

Std. Dev./ COV* 

R. 1 

LN 

>^0 * 

V <2 * 

S 

LN 

1 CM?* 

0 ^ 2 * 


















RESULTS: 



Probability 
of Failure 

Total 

CPU-Time^) 

Sample 
Size, K' c ' 

Exact (a) 

Wu/FPI 

6-W E-4 



Monte Carlo 
Variance 

Reduccion(d) 

5-072 £-f 

13-681 

1 1 589 


x 




Notes: 

(a) Exact value using POFAIL if tvo variables. If more than two, 
Wu/FPI is used; the exact should be within 5Z of this value. 

(b) CYBER 175 

(c) The number of for variance reduction and the number of for 
conventional. The values are not directly comparable. 

(d) 95 7, confidence intervals within - 10Z of p^ 

(e) Same confidence Interval as variance reduction. 









Table 3.12 Example of the Performance of a Variance Reduction Monte Carlo 
Program; EXAMPLE 9d 

DEMONSTRATING THE PERFORMANCE OF THE UA VARIANCE REDUCTION MONTE CARLO PROGRAM 
EXAMPLE H . 

PERFORMANCE FUNCTION: ^ = p— S 


Variable 

Type 

Mean/Median* 

Std. Dev./ COV* 

R 

LN. 



-S 

lk). 

l 0 • D* 

0 ' 2* 


















RESULTS: 



Probability 
of Failure 

Total 

CPU- Time C b > 

Sample 
Size, K< c > 

Exact 

Wu/FPI 

E'Y, 

x 

r><r 

Monte Carlo 

Variance 

Reduction^) 

>o2 36E-f 


iw? 

Monte Carlo 
Conventional 
(Bernoulli 
parameter) ( e ) 

>< 


1 2 4-o / ) / ?2 


Notes : 

(a) Exact value using POFAIL if two variables. If more than two, 
Wu/FPI is used; the exact should be within 5Z of this value. 

(b) CYBER 175 

(c) The number of for variance reduction and the number of Z for 
conventional. The values are not directly comparable. 

(d) 95Z confidence intervals within - 10X of p„ 

r E 

(e) Same confidence interval as variance reduction. 
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3 - 8 Comparison of Computer Costs of Variance Reduction and Conventional 
Monte Carlo 

Example 9a, b, c, and d vas designed to illustrate how computer costs 
increase as point probabilities become smaller, providing estimates at the 
same level of confidence. Figs. 3.5 and 3.6 show the relationship between 
CYBER 175 CPU execution time and the probability level for the conventional 
Bernoulli and the variance reduction estimates, respectively, for Example 9. 
Then Fig. 3.7 demonstrates how much more efficient is variance reduction 
for this problem. It should be noted that Figs. 3.5 through 3.7 relate 
onl y t0 Example 9 and cannot be presented as being characteristic of the 
relative behavior of the two methods. 

3.9 Conclusions on Variance Reduction 

Some general conclusions based on experiences exercising VARED are, 

1. Variance reduction seems to outperform conventional Monte Carlo 

\ 

consistently. However, in some esses the improvement is dramatic, in some 
cases it is modest. 

2. Related to item 1, it is difficult to predict computer costs. 

At a given confidence level, CPU time depends strongly upon the form of 
the performance function, the distribution of the variables, as well as 
the probability level. 

3. To construct a CDF, it is necessary to obtain several point proba- 
bility estimates, as it is using FPI. Thus, the variance reduction Monte 
Carlo method is not particularly effective when it is required to construct 
a distribution function of a response variable. 







BERNOULLI CPU TIME 


variance Reduction cfu.tihe 


I l I ! M If 


I I i 1 1 1 'll 


PROBABILITY OF FAILURE 

Fig. 3.7 Ratio of Bernoulli to variance reduction CPU execution time 

for CYBER 175 for point probability estimate; Example 9; 
a - 52; v - 102 
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APPENDIX A. RANDOM SAMPLES FROM GIVEN DISTRIBUTIONS 

Following are the algorithms used to generate random variates from 
the normal, lognormal, Weibull, extreme value (Type I), and the Frechet 
distributions. The computer, using a congruential algorithm, samples 
random numbers Ui from a uniform distribution U(0,1). Forms given below 
transform uniform variates to variates X ± of other models. 

Antithetic variates x^ (defined as having a negative correlation to 
x.) are generated as shown. These antithetic variates are used in the 
variance reduction method described in Section 3.0. 

A. Normal distribution, N(u, o) ; sample two uniform variates, 
and u i+r Use the Box-Muller algorithm f i, 2 J. 

X i *[> / -2 Ln(u i )’ cos(2t u 1+1 )j o + 

x! = -x, + 2u 
x i 

'V 

B. Lognormal distribution, LN(X, C^); sample two uniform variates, 
u ^ and u i+l' Use the Bo *-Muller algorithm [1, 2 ]. 

o x - /in (1 + C x ) 

U X = £n * 

x. =[exp(/-2 in(u.) cos ( 2 * u i+1 >] a x + ^ 
x' - evp(-x i + 2 p x ) 

C. Weibull distribution 

a 

Y*) ■ 1 - exp ( - (j) ) ■ u 'v l ; [0,l] 
at 

1 - u - exp (|) ) \ U[0,1] 

a 

- in a - uwf) 

x A - 8(- in (1 - u ± ) 1/a 
x| « 6(- in (u^) 1 ^ 


Thus , 
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D. EVD distribution 

F x (x) - exp (-exp(-a(X - 6)) - u % U[0,1] 
exp (-a(x - 6)) - - £n u 
- a(x - 8) - £n (- £n u) 

Thus , 

x i " B - £ £n(- £n( Ui )) 

x i “ 6 " a £n(_ £n(1 ' u i>> 

E, Frechet distribution 

F x (x) - exp - u -v, U[0,1] 

[~) m ' *n(u) 

Thus , 

x i ■ v (- £n(u.))~ 1/k 
*1 ■ v (- £n(l - Ui ))' 1/k 
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APPENDIX B. LISTING OF CONVENTIONAL MONTE CARLO PROGRAM (COMOC) 

This version runs on the VAX and the CYBER 175. It is not interactive. 

The performance function is introduced in subroutine LSFMC as XA. 

See listing. 

Card 1 Limit state function (not used in program; only printed on output) 
Card 2 Number of trials; number of variables (free format) 

Card 3 PLOT and ISTD type 

PLOT: Y^'s are sorted to construct empirical CDF 

0 « no sort 

1 » Y^'s are sorted 

ISTD; option to enter standard deviations or coefficients of 

deviations or coefficients of variation of each variable 
(if lognormal, always use COV) . 

0 « COV 

1 = Std. dev. 

Now enter each variable, its distribution type, and its moments. 

Card 4 Variable name. 

Card 5 Distribution, mean, and standard deviation 

1 = WEI (Weibull) 

2 * NORM (Normal 

3 * EVD (Extreme value distribution) 

4 * LN (Lognormal; always use median and COV) 

5 - FRE (Frechet) 


Then repeat 4 and 5 for all of the other variables. 
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1 : * PROGRAM GMC 

2s IMPLICIT DOUBLE PRECISION <A-H,0-Z> 

3s DIMENSION I NAME (40) f XMEAN (40) , XSTD (40) , DIST (40) , DTRANS (40) , X 

4s DIMENSION Y (30000) ,F (5) ,AL (40) , BE (40) 

5s COMMON /TWO/ PI f PI2 


6s 
7s 
8s 
9s 
10s 
11s 
12: 
13s 
14s 
15s 
16s 
17s 
18s 
19s 
20 s 
21s 

22 s 

23 s 
24s 

25 s 

26 s 

27 s 

28 s 
29: 
30: 
31s 
32: 
33s 
34; 
35s 
36 s 
37s 

38 s 

39 s 
40: 
41s 
42 s 
43s 
44 s 
45: 


CHARACTER*80 BRS , F I N , FOUT , DTRANS*7 , AA*7 , BB«6 , CC*3 , DD#3 , EE* 
DATA AA/ ‘WE I BULL'/ 

DATA BB/ 'NORMAL'/ 

DATA CC/'EVD'/ 

DATA DD/'LOB’/ 

DATA EE/ 'FRECHET'/ 

651 FORMAT (A10) 

8004 CONTINUE 

READ <5 , ' (A) ' ,END*8B8B> BRS 
READ (5,*) K,N 

READ (5,*) PLOT 1 ,PL0T2 , 1 STD 
READ (5,*) I SEED 
DO 7901 I-1,N 
READ (5,651 ) I NAME ( I ) 

READ (5,*) DIST (I) , XMEAN ( I ) v XSTD ( I > 

7901 CONTINUE 

IF ( ISTD. EQ. 0) THEN 
DO 913 I*1,N 

IF (DIST ( I ) . EQ. 4. ) BO TO 913 
XSTD ( I ) *XMEAN ( I > *XSTD ( I ) 

913 CONTINUE 

END IF 
C 
C 

DO 1234 1*1 ,N 
AL ( I ) *0. D0 
BE ( I ) *0. D0 
1234 CONTINUE 

PI *4. D0*DATAN ( 1 . D0) 

PI2-PI+PI 
DO 1 I * l , N 
IF (DIST < I ) .EQ. 1.) DTRANS ( I ) *AA 
IF (DIST (I) .EQ.2. ) DTRANS ( I >*BB 
IF (DIST ( I ) . EQ. 3. ) DTRANS ( I > *CC 
IF (DIST ( I ) . EQ. 4. ) DTRANS ( I >*DD 
IF (DIST (I) . EQ. 5. ) DTRANS ( I > *EE 

IF (DIST (I) .EQ. 1. > CALL WEI (XMEAN ( I > , XSTD ( I > , AL ( I ) ,BE ( I ) > 

IF (DIST ( I ) . EQ. 3. ) CALL EVD(XMEAN(I> , XSTD(I) ,AL(I) , BE(I) ,PI) 
IF (DIST (I) . EQ. 5. ) CALL FRE (XMEAN ( I ) , XSTD ( I > , AL ( I ) ,BE ( I ) ) 

1 CONTINUE 


CONVENTIONAL 
MONTE CARLO 
PROGRAM (COMOC): 
Runs on the VAX 
or CYBER 175 


46: C 

47: C* THE DATA IS PRINTED OUT. 


48s C 

49 s 

50 s 
51s 
52: C 
53 s 
54: 

55 s 

56 s 

57s 3 

58: 

59 s 


WRITE (6 f 1 1 ) 6RS,K,N 
WRITE (6, 12) 

WRITE (6, 13) ( I NAME ( I > , DTRANS (I) ,XMEAN(I) ,XSTD(I> , I-1,N> 
GENERATE RANDOM # AND CORRESPONDING RANDOM VARIABLE 
NUM*0 

DO 4 1*1, K 
DO 3 J* 1 , N 

CALL GENX (DIST(J) ,AL(J) ,BE(J) ,X(J) ,XMEAN(J) ,XSTD(J) ,ISEED) 
CONTINUE 

CALL LSFMC < Y ( I ) ,N,X> 

IF(Y(I) .LE.0.D0) NUM-NUM+1 
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60: 
61 : 
62: 
63: 
64: 
65: 
66 : 
67: 
68 : 
69: 
70: 


12 

12 


C 

c* 

c* 

c 


71: 

C 

72: 

C* 

73: 


74: 


75: 

92 

76: 


77: 

1017 

78: 


79: 


80: 

3030 

81: 

1003 

82: 


83: 


84: 


85: 


86: 


87: 


88: 


89: 


90: 

3031 

91: 


92: 

67 

93: 


94: 


95: 


96: 


97: 


98: 


99: 


100: 


101: 


102: 


103: 

1009 

104: 

3456 

105: 


106: 

1 1 

107: 


108: 


109: 

12 

1 10: 


111: 

13 

112: 

15 

113: 


114: 

- 

115: 

17 

116: 

19 

117: 


118: 1 

B88S 

119: 

125 


4 CONTINUE 

3 CALL STAT (Y ,K, YMEAN, Y5TD, YMED, YCOV) 

4 WR I TE (6 , 15) YMEAN ,YSTD, YMED, YCOV 

ROUTINE TO ACCUMULATE NUMBER OF TRIALS WITH NEGATIVE Y(I> 
VALUES AND PRINT OUT RESULTS Vm 

RATIO « DBLE(NUM) /DBLE(K> 

WRITE (6,9) NUM, RATIO 
9 FORMAT ( / , 10X , 'NUMBER OF NEG 
+' PERCENT OF TRIALS* ' ,F9. 6) 


Y VALUES * ' 1 15, * . ' ,4X , 


THE SORTED VALUE OF Y AND THE EMPIRICAL CDF ARE PRINTED 
I F < PLOT 1 . EO . 0 . ) GO TO 92 
CALL OSORTCY.K) 

IF (PL0T2. EO. 0. ) GO TO 3456 
WRITE<6, 1017) 

FORMAT (////, 14X, 'SORTED VALUES OF Y AND THE EMPIRICAL CDF 
J 1 * 1 

J2*5 

WRITE (6, 1003) Jl, (Y(I) ,I-J1, J2> 

FORMAT (IX, 'I - ' , 15 ,5E13. 5) 

J 1*J 1+5 
J2-J2+5 

IF(Jl.GT.K) GO TO 3031 
IF ( J2. GT . K) THEN 
J2=K 

GO TO 3030 
END IF 
GO TO 3030 
CONTINUE 
WR I TE <6 , 67 ) 

FORMAT </) 

J=0 
Jl = l 

DO 1009 1=1, K 
J=J + 1 

F ( J) = (DBLE ( I ) -. 5) /DBLE (K) 

IF < J. EO. 5. OR. I . EQ. K) THEN 
WRITE (6, 1003) Jl , <F(L) ,L*1,J) 

J=0 

J 1 =J 1 +5 
END IF 
CONTINUE 
CONTINUE 


LIMIT STATE FUNCTION : ' , A ,5 C / ) , 10X , 

+ SAMPLE SIZE, K* ', I7//10X, 'NUMBER OF RANDOM VARIABLES, N=',I3/, 
FORMAT <26X , 'RANDOM VARI ABLES //10X VARIABLE 2X , 

+ 'DISTRIBUTION ' ,8X, 'MEAN', 12X, 'STD DEV') 


+STD DEV *' ,E13. 5//10X, 'MEDIAN * ' ,E13. 5// 10X , 'COV 


,E13. S//10X 


GO TO 8004 
CONTINUE 



£nd 

SUBROUT I NE ST AT < U , M , XM , STD , XMED , COV ) 
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120; 

121i 

122; C 
123: C* 
124: C* 
125: C 
126: 
127: 
128: 
129: 
130: 
131: 
132: 
133: 
134: 
135: 
136: 
137: 
138: 
139: 
140: 
141: 
142: 

143: 

144: 

145: 

146: 

147: C 
148: 

149: 

150: 

151: 1 

152: 
153: 2 
154: 
155: 
156: 
157: 3 
158: 
159: 4 
160: 
161: 
162: 

163: 

164: 

165: 5 
166: 

167: 

168: 

169: 

170: C 

171: C 

172: 

173: 

174: 

175: 
176: 7 
177: 
17B: 

179: 


THIS SUBROUTINE CALCULATES THE STATISTICS (MEAN, STD DEV, MED I AN 
OF Y FUNCTION. 


IMPLICIT DOUBLE PRECISION (A-H t O-Z) 

DIMENSION U (M) 

XK-M 

XM-0. 

DO 63 I-1,M 
XM»XM+U<I) 

63 CONTINUE 
XM-XM/XK 
STD-0. 

* DO 64 I-1,M 

STD-STD* (U ( I ) -XM) **2 

64 CONTINUE 
STD-STD/ ( XK- 1 . D0 > 

STD-DSQRT (STD) 

COV-STD/XM 

XMED-XM/DSQRT ( 1 . D0+COV**2) 

RETURN 

END 

SUBROUTINE GENX (D 1ST, ALPHA, BETA, X,XMEAN,XSTD, I SEED ) 
IMPLICIT DOUBLE PRECISION (A-H,0-Z> 

COMMON /TWO/ PI ,PI2 

IDIST-INT (DIST+. 1) 

AA-RAN ( I SEED > 

GO TO <1,2, 3, 4, 5), IDIST 
X-BETA* (-DLOG <AA> )**<!. D0/ALPHA) . 


GENX obtains 
random samples 
from distributions 


RETURN 

BB-RAN ( I SEED) 

E-D5QRT (-2. D0*DLOG (AA) > 

X«E*DCOS <PI2*BB) *XSTD+XMEAN 
RETURN 

X-BETA-DLOG (-DLOG (AA) ) /ALPHA 
RETURN 

BB-RAN < I SEED) 

SDX-DSQRT (DLOG < 1 . D0+X5TD**2) > 
UX-DLOG ( XMEAN) 

E-DSQRT ( -2 . D0*DLOG ( AA ) ) 

X-DEXP (E*DCOS (PI2*BB) *SDX+UX) 
RETURN 

X -BET A* < -DLOG (AA) ) ** (-1 . D0/ ALPHA) 


RAN is library 
uniform random 
number generator 
for CYBER 175 


RETURN 


£feU2 

SUBROUT I NE B I SECT ( COV , I S I GN , ALPHA ) 

IMPLICIT DOUBLE PRECISION <A-H,0-Z> 

ISIGN - If WE I BULL DIST. 

■ 2) FRECHET DI8T. 

F (X ,COV) ■- ( 1 . D0+COV**2) *GAMMA (X) **2+GAMMA <2. *X ) 
IF (ISIGN. EQ. 1 ) X 1 «COV** ( 1 .08) 

I F < I S I GN . EO . 2 ) XI -COV** ( . 677 ) /2 . 33 

IF ( I S I GN . EQ . 2 . AND .XI. GT . . 49D0 ) X1-. 48999999 

IF ( ISIGN. EQ. 1 ) F 1 »F < X 1 , COV ) 

IF ( ISIGN. EQ. 2) Fl-F < — X 1 ,COV) 

IF (DABS (FI ) . LE. 1 . D-10) GO TO 1 
X2-X1+.01D0 
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180: 

181 s 
182: 

183: 

184: 

185: 

186: 

187: 20 
188: 2 
189: 

190: 

191: 

192: 

193: 

194: 

195: 1 
196: 

197: 

198: 

199: 

200 : 

201 : 

202 : 

203: 

204: 

205: 

“MSi 

207: 

208: 

209: 

210 : 

211 : 

212 : 

_213j 

214: 

215: 

216: 

217: 

218: 

219: 

220 :' 

221 : 

222 : 

223: 

224: 

225: 

226: 

227: 

228: 

229: 

230: 

231: 

232: 

233: 

234: 

235: 

236: 

237: 

238: 

239: 


456 


* 

* 


3 

457 


IF ( ISIGN. EQ. 1 ) F2-F (X2.C0V) 

IF ( ISIGN. EQ. 2) F2«F(-X2,C0V> 

F12«F1*F2 

IF <F12. LT . 0. ) SO TO 20 
IF (DABS (FI ) . GT. DABS <F2) ) X1-X2 
IF (DABS <F1 ) .LT. DABS (F2) > X1-X1-.01DO 
GO TO 7 
CONTINUE 
X3«(X1+X2)*.5D0 

IF ( ISIGN. EQ. 1 ) F13-F (XI ,COV> #F (X3.C0V) 

IF ( ISIGN. EQ.2> F13«F<-X1,C0V)*F<-X3,C0V> 

IF (F13. LT .0. ) X2-X3 
IF (FI 3. GT . 0. ) X1-X3 
DX«DABS(X1-X2> 

IF(DX.GE. l.D-9) GO TO 2 
ALPHA- 1. DO/ XI 
RETURN 

END 

SUBROUTINE WEI (XMEAN, XDEV, ALPHA, BETA) 

IMPLICIT DOUBLE PRECISION (A-H f O-Z> 

COV-XDEV/ XMEAN 
CALL B I SECT ( COV , 1 , ALPHA ) 

AL 1-1. DO/ ALPHA 
BETA-XMEAN/GAMMA < ALl ) 

RETURN 
END 

"SUBROUTINE FRE (XMEAN, XDEV .ALPHA, BETA) 

IMPLICIT DOUBLE PRECISION <A-H,0-Z> 

COV-XDEV/ XMEAN 

CALL BISECT (COV, 2 .ALPHA) " \ 

AL 1-1. DO /ALPHA 
BETA= XMEAN /GAMMA ( -AL 1 ) 

RETURN 

END 

SUBROUTINE EVD < XMEAN , STD , ALPHA , BETA , PI ) 

IMPLICIT DOUBLE PRECISION <A-H,0-Z) 

ALPHA-PI/ <STD*DSQRT (6. DO) ) 

BETA-XMEAN-. 5772 1566490 153/ ALPHA 
RETURN 
END 

DOUBLE PRECISION FUNCTION GAMMA (Yl) 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

COMMON /TWO/ PI, PI 2 
X-Y1+1.D+0 
Z«X 

IF (X. GE. 6. 0D-t-0) GO TO 456 
N*INT(X) 

Z* (6. 0D+0) -N+X 
Y-l.D+0/Z**2 

ALG- ( Z-. 5D+0) *DLOG ( Z ) ♦. 5D+0*DL0G (PI2) - 
Z-d.D+0/ (12.D+0*Z) )*( ( (Y/0. 14D+3-1 . D+0/0. 105D+3)*Y+ 
1 . D+0/. 3D+2) *Y-1 . D+0) 

IF(X.GE.6.D+0)G0 TO 457 
ITE-6-N 
DO 3 J«t , ITE 
A-X+J-i.D+0 
ALG-ALG-DLOG (A) 

CONTINUE 

GAMMA-DEXP(ALG) 

RETURN 


The gamma 
function 


Computes 

EVD 

parameters 


Computes 

Frechet 

parameters 


Computes 

Weibull 

parameters 


BISECT used to 
compute Weibull 
and Frechet 
shape parameter 
(exponent) 
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273 

274 

275 

276 

277 

278 

279 

280 
281 
282 

283 

284 

285 

286 

287 

288 
289 1 


END 

SUBROUTINE QSORT<A f N> 

IMPLICIT DOUBLE PRECISION <A-H t O-Z> 
DIMENSION A(N) ,KSL (24) V KSR(24) 

KS-l 

KSL <1 ) • 1 
KSR ( 1 ) »N 
10 CONTINUE 

L-KSL(KS) 

KR-KSR(KS) 

KS-KS-1 

20 CONTINUE 

I-L 
J-KR 

LR-<L+KR>/2 
X-A(LR) 

30 CONTINUE 

IF ( A ( I ) . LT . X ) THEN 
I-I + l 
60 TO 30 
END IF 

40 CONTINUE 

IF < X . LT . A ( J > ) THEN 
J-J-l 
GO TO 40 
END IF 

IF(I.LE.J) THEN 
W«A ( I ) 

A < I ) ■ A ( J ) 

A ( J ) «W 
I*I + 1 
J=J-1 
END IF 

IF(I.LE.J) 60 TO 30 
IF(I.LT.KR) THEN 
KS=KS+1 
KSL (KS>«I 
KSR <KS) “KR 
END IF 
KR=J 

IF (L. LT. KR) 60 TO 20 
IF <KS. NE. 0) 60 TO 10 
RETURN 

END 

SUBROUTINE LSFMC < XA,N, X > 

IMPLICIT DOUBLE PRECISION <A-H t O-Z 
DIMENSION X <N> 

XA»X(1)-X(2) ■<] 

RETURN 
END 


This is where the 
limit state is 
introduced 


This is the 
sort routine, 
QUICKSORT 
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APPENDIX C. THE SORT ROUTINE: "QUICKSORT" 

QUICKSORT is described in detail in the book by Wirth [7], who describes 
its performance as "spectacular," and claims that it is the best sorting 
method on arrays known so far. The method is based on exchanges and the 
inventor C.A.R.Hoare recognized that sorting becomes most efficient when 
exchanges are made over large distances. 

The table below shows execution times (in milliseconds) consumed by 
several proposed sorting methods as executed by the PASCAL system on a 
CDC 6400 computer. The three columns contain times used to sort the 
already ordered array, a random permutation, and the inversely ordered 
array. The left figure in each column is for 256 items, and right one 
for 512 items. 

In summary, the computational effort needed for QUICKSORT is of the 
order of n log n. 


Ordered Random Inversely Ordered 


Straight insertion 

Binary insertion 

Straight selection 

Bubblesort 

Bubblesort with flag 

Shakersort 

Shellsort 

Heapsort 

Quicksort 

Mergesort 


12 

23 

366 

56 

125 

373 

489 

1907 

509 

540 

2165 

1026 

5 

8 

1104 

5 

9 

961 

58 

116 

127 

116 

253 

110 

31 

69 

60 

99 

234 

102 


1444 

704 

2836 

1327 

662 

2490 

1956 

695 

2675 

4054 

1492 

5931 

4270 

1645 

6542 

3642 

1619 

6520 

349 

157 

492 

241 

104 

226 

146 

37 

79 

242 

99 

232 


Execution Times of Sort Programs. 
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APPENDIX D. LISTING OF THE VARIANCE REDUCTION MONTE CARLO PROGRAM (VARED) 

This version runs on the VAX and the CYBER 175. It is not interactive. 

The performance function is introduced in subroutine LSFMC, then compiled 
and linked to the rest of the program. 



Any integer number between 0 and 262,139 to start the random sampling. 
Ex: 23, 579, etc. 

Enter variable name. (Free format) 


Card 5 


E- 59 


Card 6 Enter corresponding distribution, mean, and standard deviation 
(if LN always input median and COV); Ex: 1, 20, 2 

a. dist . * 1 * Weibull 

2 ■ Normal 

3 - EVD 

4 ■ Lognormal (LN) 

5 * Frechet 


Then repeat 5 and 6 for all of the other variables. 
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1 : 


4: 
5 : 
6 : 
7: 
8 ; 
9: 
IB): 
11 : 
12 : 
13: 
14: 
15: 
16: 
17: 
18: 
19: 
28: 
21 : 


24 

25 

26 

27 

28 
29 
38 
31 


34 

35 

37 

38 

39 

48 

41 

42 

43 

44 

45 

46 

47 

48 

49 
58 

51 

52 

53 

54 

55 

56 
57: 
56: 
59: 


* PROGRAM GMC 

IMPLICIT DOUBLE PRECISION <A-H,0-Z) 

DIMENSION I NAME (28) ,XMEAN(2B> ,XSTD<28> ,DIST(2B) ,DTRANS (28) ,XC 
DIMENSION Y(1008B> ,F<5) ,AL(28> ,BE<28> ,XA<28> ,TXC2B> ,TS<28> 
COMMON /TWO/PI, SPI2,PI2 

CHARACTER*70 GRS , F I N , FOUT , AA*7 , BB*6 , CC*3 , DD*3 , EE*7 
CHARACTER*7 INAME, DTRANS 
DATA AA/ 'WE I BULL’/ 

DATA BB/' NORMAL'/ 

DATA CC/ 'EVD '/ 

DATA DD/ 'LOG'/ 

DATA EE/ FRECHET '/ 

C 

READ (5, ' (A) ' ,END=BB88) GRS 
READ (5,*) K , N , EPS 
READ ( 5 , •* ) ZAL, GAM, I STD, PLOT 
READ (5,*) I SEED 
DO 7981 I = 1,N 
READ (5, ' (A) ' > I NAME ( I > 

READ (5,*) DIST(I) , XMEAN ( I ) ,XSTD<I> 

7981 CONTINUE 

8884 CONTINUE 

IF(ISTD.EQ.B) THEN 
DO 913 1=1, N 

IF (DIST ( I ) . EO. 4. ) GO TO 913 
XSTD < I ) =XMEAN ( I ) *X5TD ( I ) 

913 CONTINUE 

END IF 

IF (K.GT. 18888) K=18888 
C 
C 

DO 1234 I = 1 , N 
AL ( I ) =8. D8 
BE ( I > =8. D8 

IF (DIST < I ) . EQ. 4. ) THEN 

TX ( I ) = XMEAN ( I ) *DSQRT ( 1 . DB+XSTD ( I ) **2> 

TS ( I ) =TX < I ) *XSTD( I ) 

ELSE 

TX ( I ) =XMEAN < I ) 

TS ( I > =XSTD ( I ) 

END IF 

1234 CONTINUE 

P I = 4 . D8*DATAN ( 1 . DB) 

P I 2=P I +F I 

SF 12=1. D8/DS0RT (PI2) 

DO 1 I = 1 , N 

IF (DIST ( I ) . EQ. 1 . > DTRANS ( I ) =AA 
IF (DIST ( I ) . EQ. 2. ) DTRANS ( I > =BB 
IF (DIST < I ) . EQ. 3. ) DTRANS ( I > =CC 
IF (DIST (I) . EQ. 4. ) DTRANS ( I ) =DD 
IF (DIST ( I ) . EQ. 5. ) DTRANS ( I > *EE 

IF(DIST(I> .EQ. 1. ) CALL WEI ( XMEAN ( I ) , XSTD ( I ) , AL ( I ) , BE ( I ) ) 

IF (DIST ( I ) . EQ. 3. > CALL EVD ( XMEAN ( I ) , XSTD ( I ) , AL < I ) , BE ( I ) , PI ) 
IF (DIST ( I ) . EQ. 5. ) CALL FRE ( XMEAN ( I ) , XSTD ( I ) , AL ( I ) , BE < I ) ) 

1 CONTINUE 
C 

C* THE DATA IS PRINTED OUT. 

C 

C MAIN LOOP USING ANTITHETIC VARIANCE REDUCTION METHOD 


Program VARED. Monte 
Carlo using variance 
reduction method; runs 
on the VAX or CYBER 175 




62 : 
61: 
62: 
63: 
64: 
65 : 
66 : 
67: 
68 : 
69: 
70: 
71: 
72: 
73: 
74: 
75: 
76: 
77: 
78: 
79: 
80: 
81: 
82: 
8 j: 

84: 

85: 

86 : 

87: 

88 : 

89: 

90: 

91: 

92: 

93: 

94: 

95: 

96: 

97: 

98: 

99: 

100 : 

101 : 

102 : 

103: 

104: 

105: 

106: 

107: 

108: 

109: 

110 : 

111 : 

112 : 

113: 

114: 

115: 

116: 

117: 

118: 

119: 


701 

700 

C 


96 

559 

561 

563 


98 


70; 


FIND MAX. IMPACT VARIABLE 
DG=0.D0 

CALL LSFMC (G,N,TX) 

DO 700 1=1, N 
TX(I)«TX(I)+TS(I) 

CALL LSFMC <DGB,N,TX> 

DGA=DGB-G 

IF (DABS (DGA) .LE. DABS (DG) ) GO TO 701 

IV=I 

DG=DGA 

TX ( I ) =TX ( I ) -TS ( I > 

CONTINUE 
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' ,F6. 2, ' X',//, 


X ( * , 12, ') ',/) 


702 

C 

12: 


WRITE (6,11) GRS,K,N 
WRITE (6,96) ZAL ,GAM 
FORMAT ( 10X , ' CONFIDENCE INTERVAL 
r 10X, 'GAMMA = ' , F6. 2, // /) 

WRITE (6,559) IV 

FORMAT ( 10X , 'MAX. IMPACT VARIABLE ■ 

IF (DG.LE.0.D0) WRITE (6,561 > 

FORMAT (10X, 'VARIABLE TYPE IS STRESS',///) 

IF (DG.GT.0.D0) WRITE (6,563) 

FORMAT ( 10X , 'VARIABLE TYPE IS STRENGTH',///) 

WRITE (6,12) 

WRITE <6, 13) ( I NAME ( I ) , DTRANS ( I ) , XMEAN ( I ) , XSTD ( I > ,1 = 1 ,N) 

CALCULATE PROB. OF FAILURE 

K 1 = 1 

K2=K 

ICO=l 

CONTINUE 

DO 702 I=K 1 , K2 

DO 703 J= 1 , N 

IF ( J . EO. IV) GO TO 703 

CALL GENX (DIST ( J ) , AL ( J > ,BE ( J ) , X ( J > , XA ( J ) , XMEAN (J) , XSTD( J> , ISE 
CONTINUE 

IF (DG. GT. 0. D0) A=TX ( I V) -3. D0*TS ( IV) 

IF (DG. LE. 0. D0) A=TX ( I V) +2. D0*TS (IV) 

B=A-*-TS ( IV) 

CALL SECA < EPS ,A,B, I V , N , X ) 

CALL CDFPDF (DIST ( IV) , AL ( I V) ,BE ( IV) , X ( I V) , XMEAN (IV) ,XSTD(IV> , 

* 1, CDF 1, PDF) 

IF (DG. LE. 0. D0) CDF 1=1. D0-CDF 1 
IF (DG. GT.0. D0) A=TX ( IV) —3. D0*TS (IV) 

IF (DG. LE. 0. D0> A=TX ( IV) +2. D0*TS (IV) 

B=A+TS (IV) 

CALL SECA (EPS ,A,B,IV,N,XA> 

CALL CDFPDF (DIST ( IV) ,AL(IV) ,BE(IV) ,XA(IV> , XMEAN (IV) , XSTD (IV) , 

* 1 , CDF2 , PDF ) 

IF (DG. LE. 0. D0) CDF2=1 . D0-CDF2 
Y < I ) = ( CDF 1 +CDF2 > * . 5D0 
CONTINUE 

CALL STAT (Y,K1 , K2 , YMEAN , YSTD , YMED , YCOV > 

IF(ICO.EO.l) THEN 

YM=YMEAN 

YS=YSTD 

YME=YMED 

YC=YCOV 

YM 1 =YM 

ELSE 
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120 : 
121 : 
122 : 
123: 
124: 
125: 
126: 
127: 
128: 
129: 
130: 
131: 
132: 
133: 
134: 
135: 
136: 
137: 
138: 
139: 
140: 
141: 
142: 
143: 
144: 
145: 
146: 
147: 
148: 
149: 
150: 
151 : 
152: 
153: 
154: 
155: 
156: 
157: 
158: 
159: 
160: 
161 : 
162: 
163: 
164: 
165: 
166: 
167: 
168: 
169: 
170: 
171 : 
172: 
173: 
174: 
175: 
176: 
177: 
178: 
179: 


< K*VM 1 + (K2-K) * YMEAN) /K2 

V51 =y 5**2* (K-l > +K*YM1**2+YSTD**2* (K2-K-1 ) + (K2-K) *YMEAN**2 
YS2= YS 1 — K2*YM**2 
YS=D5QRT ( YS2/ (K2-1 ) ) 

YC=YS/YM 

YME=YM/DSQRT ( 1 . D0+YC**2 ) 

END IF 

ZAL1=. 005D0* (100. D0+ZAL) 

Z AX=X INV ( ZAL 1 ) 

ZX*ZAX*YC/DSDRT (D&LE (K2) ) 

PL=YM*(1.D0-ZX) 

PU=YM* ( 1 . D0+ZX) 

WRITE (6, 176) YM,ZAL,PL,PU 

176 FORMAT (///, 10X, 'ESTIMATE OF P - ',E13.5,//, 

# 10X,F5.2,' X CONFIDENCE INTERVALS ARE',//, 

* 10X,'PL = ' ,E13.5,5X, 'PU « ' ,E13. 5, ///) 


C 

WRITE (6, 15) YMEAN, YSTD , YMED, YCOV 
IF (PLOT. EQ. 0. ) GO TO 3456 



J 1 = 1 
J2=5 

3030 WRITE (6, 1003) Jl, (Y(I) ,I=J1,J2> 

1003 FORMAT (IX, 'I = ’,I5,5E13.5) 

J 1 =J 1 +5 
J2=J2+5 

I F ( J 1 . GT . K2 ) GO TO 3031 
IF ( J2. GT. K2> THEN 
J2=K2 

GO TO 3030 
END IF 
GO TO 3030 

3031 CONTINUE 

WR I TE (6 , 67 ) 

67 FORMAT (/) 

J=0 
J 1 = 1 

DO 1009 1=1, K2 
J=J+1 

F ( J ) = ( D&LE ( I ) -. 5) /DBLE (K2) 
IF(J.E0.5.0R. I.EQ.K2) THEN 
WRITE (6, 1003) Jl , (F (L) ,L=1 , J) 
J=0 


J 1 = J 1 +5 
END IF 

1009 CONTINUE 

3456 CONTINUE 

Kl=k>l 

K2= (YC*ZAX/GAM) **2+1 
IF(ICO.EQ.l) WRITE (6,99) GAM , K2 
99 FORMAT (//, 10X , 'K FOR GAMMA = ',F6.2,' IS ',16) 


ICO=ICO+l 

IF ( ICO. EQ. 2. AND. K2. GT . K> GO TO 98 
11 FORMAT ( 1H1 ,5 (/) ,30X , 'MONTE CARLO SOLUT ION , 5 ( / ) , 1 0X , 
+ 'LIMIT STATE FUNCTION : ' , A,5 ( / ) , 10X , 

+' SAMPLE SIZE =', 17// 10X , 'NUMBER OF RANDOM VARIABLES 


13 //) 



630 1 


FORMAT (26X , ‘RANDOM VARIABLES ' ,//10X , ’VARIABLE ,2X , 
S^STRIBUlioN- ,SX, ’ MEAN » 12X , ' STD DEV ') 

13 F0RMAT</llX,A7 f 5X,A7,5X,El2.5,5X,E12.5> 

15 FORMAT (Z////10X, 'STATISTICS OF P I //W* MEAN - 
+'STD DEV *’ ,E13.5//10X, 'MEDIAN » f E13.5,//10X, COV 

+E1 IF^ANSl .EQ. 'F ' .OR. ANS1.EQ. '♦ * > BO TO 8300 

FORMAT < ’ DO YOU HAVE ANOTHER DATA SET ?<Y/N> 

READ (5,8001) ANS3 __ __ 

IF(ANS3.E0. ' Y ' . OR. ANS3. EQ. y ) BO TO 8304 

CONTINUE 
STOP 
END 
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,E13.5//10X, 




SUBROUTINE SECA (EPS , A ,B , I V ,N , X ) 
IMPLICIT DOUBLE PRECISION (A-H,0-Z> 
DIMENSION X (N) 
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240: 

241: 

242: 

243: 

244: 1 

245: 

246: 

247: 

246: 

249: 

250: 

251 : 
252: 
253: 



255: 
256: 
257: 
258: 
259: 
260: 
261: 
262: 
263: 
264: 
265: 
266: 
267: 
26B: 
269: 
270: 
271: 
272: 
273: 
274: 
275: 
276: 
277: 
278: 
279: 
280: 
281: 
282: 
283: 
284: 
285: 
286: 
287: 
288: 
289: 
290: 
291: 
292: 
293: 5 
294: 
295: 
296: 
297: 

. 

299: 


X ( IV) =A 

CALL LSFMC (U , N , X ) 

X ( IV) =B - 1 

CALL LSFMC (V,N,X) 

CONTINUE 

IF (DABS (X (IV) -A) .GE. EFS) THEN 
X(IV)*B-V*<B-A)/(V-U> 

A=B 

B=X (IV) 
u=v 

CALL LSFMC (V,N, X) 

GO TO l 
END IF 
RETURN 
END 


This defined the performance function 


This subroutine determines the 
point at which the CDF is 
evaluated for the maximum 
impact variable 


SUBROUT I NE CDFPDF ( D I ST , ALPHA , BETA , X , XMEAN , XDEV , I CDF , CDF , PDF ) 
IMPLICIT DOUBLE PRECISION <A-H f O-Z> 

COMMON /TWO/PI ,SPI2,PI2 
IDIST=INT (DIST+. 1 ) 

GO TO ( 1 ,2, 3, 4 , 5 ) , IDIST 
1 RB=X /BETA 


EW=RB** ALPHA 

IF <EW. GT. 200. ) EW=200 . 

EXF'WE I =DEXP ( -EW ) 

CDF= 1 . D0-EXPWE I 

IFdCDF.EQ. 1) GO TO 10 

PDF= (ALPHA/BETA) * (EW/RB) *EXPWEI 


Evaluates 
the CDF 


GO TO 10 

2 Z= (X-XMEAN) /XDEV ( 

CDF=CDFNOR (Z) 

IF ( ICDF. EQ. 1 ) GO TO 10 
PDF=SP I 2*DEXP ( -Z**2* . 5D0) /XDEV 


GO TO 10 

3 EE=ALPHA* ( x-BETA) 

IF (EE. GT. 200. ) EE=200 . 

YV=DEXP (-EE) 

IF ( YY . GT . 200 . ) YY=200. 

CDF = DEXF' (-YY) 

IF ( ICDF. EQ. 1 ) GO TO 10 
EY=EE+ YY 

IF <EY. GT. 200. ) EY=200. 
PDF=ALPHA*DEXP (-EY) 

GO TO 10 

4 CX2l=XDEV**2+l . D0 
YMEAN=DLOG (XMEAN) 

YDEV=DSQRT (DL0G(CX21 ) ) 

Z= (DLOG (X) -YMEAN) / YDEV 
CDF=CDFNOR (Z) 

IFdCDF.EQ. 1) GO TO 10 
EZ=- (Z**2) *. 5D0 

IF (EZ. LE. -200. ) EZ=-200. 
PDF=SPI2*DEXP(EZ> / (YDEV*X) 

GO TO 10 

TEMP= (BETA/X ) **ALPHA 
CDF=DEXP (-TEMP) 

IF (ICDF. EQ. 1 ) GO TO 10 
PDF=CDF*TEMF* ALPHA/ X 
10 RETURN 

END — — “ 

DOUBLE PRECISION FUNCTION XINV (Z) 
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354: 
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IMPLICIT DOUBLE PRECISION (A-H f O-Z> 
F ( X ,P 1 > *P1 -CDFNOR < X ) 

Y=Z 

IF(Z.GT.0.5D0> Y=1.D0-Z 
C0=2. 515517D0 
C1*0. 802853D0 
C2*0. 010328D0 
D1-1.4327BBD0 
D2«0. 189269D0 
D3*0. 001308D0 

T- ( -2 . D0*DLOG ( Y > > ** . 5D0 
DNUM»C0+T* (C1+T*C2) 
DNOM«1.0D0+T*(Dl-t-T*(D2+T*D3) ) 

X=T- (DNUM/DNOM) 

IF ( Z . LT. 0. 5D0) X«-X 
A=X 

B*X+. 001D0 
V*F (B,Z> 

U=F (A, Z) 

XX*B 

CONTINUE 

IF (DABS (XX-A> . GE. 1 . D-10) THEN 
XX=B-V* (B-A) / (V-U) 

A=B 
B=XX 


The inverse normal 
using the secant 
method 


V=F(XX,Z> 
GO TO 1 
END IF 
XINV=XX 
RETURN 
END 


DOUBLE PRECISION FUNCTION CDFNOR(Z) 

C THIS FUNCTION COMPUTES THE NORMAL CDF. 

IMPLICIT DOUBLE PRECISION (A-H,0-Z> 

COMMON /TWO/PI ,SPI2,PI2 

DATA A/0. 31 93B153D0/,B/-0 .-3565637B2D0/ ,C/1 .7814779^700/ , 
+D/-1 . 82125597BD0/ ,E/ 1 . 330274429D0/ 

EZ=- ( Z**2> *. 5D0 
CDFNOR=0 . 0D0 

I F ( E Z . LE . -200 . 0D0) GO TO 1 
ZX=SPI2*DEXP (EZ> 

IF (DABS < Z ) . GT . 6. D0) GO TO 2 
T*1 . D0/ ( 1 . D0-*- (0. 2316419D0*DABS (Z> > ) 

CDFNOR=ZX*T* ( A+T* (B+T* (C+T* (D+T*E> > > > 

GO TO 1 

2 CDFNOR*ZX** 1 ! D8-Z2. < 1 . D0-3. D0.Z2. < 1 . D0-5.D0.Z2I ) > /DABS < Z) 

1 IF ( Z . GT. 0. 0D0 ) CDFNOR-1.0P0-CDFNOR 
, RETURN 

SUBROUTINE STAT (U,K1 , K 2 , XM ,feTD, XM&D ,Cdv> 


C* THIS SUBROUTINE 
C* OF Y FUNCTION. 


CALCULATES THE STATISTICS (MEAN, STD DEV,MEDIAN,COV> 


IMPLICIT DOUBLE PRECISION (A-H,0-Z> 
DIMENSION U(K2> 

XK=K2-K1+1 

XM*0. 
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403 

404 
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408 
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410 

411 

412 

413 

414 

415 

416 
417; 
418 
419: 


DC 63 I — K 1 , K2 

Xtf=Xh-rU ( I ) 

63 CON T I NUE 
X M * X M / X k 
STD-0 . 

DO 64 I *K 1 , K2 
STD=STD+ <U < I ) -XM) **2 

64 CONTINUE 
STD*5TD/(XK-1.D0> 

STD=DSORT <STD) 

COV=STD/XM 

XHED=XM/DSQRT ( 1 . D0+COV**2> 

RETURN 

END 

SUBROUTINE GENX (DIST, ALPHA , BET A , X , X A , XME AN , XSTD , I SEED > 
IMPLICIT DOUBLE PRECISION (A-H,0-Z> 

COMMON /TWO/PI ,SPI2,PI2 
C 

IDIST=INT (DIST+. 1 ) 

AA*RAN ( I SEED) 

GO TO (1,2, 3, 4,5), IDIST 

1 X=BETA* (-DLOG (AA) ) *# ( 1 . D0/ALPHA) 

XA=BETA* (-DLOG < 1 . D0-AA) ) #* ( 1 . D0/ALPHA) 

RETURN 

2 BB=RAN ( I SEED) 

E=DSQRT (-2. D0*DLOG (AA) ) 

X=E*DCOS <PI2*BB) *XSTD+XMEAN 
XA=-X+2. D0*XMEAN 
RETURN 

3 X=BETA-DLOG (-DLOG <AA) ) /ALPHA 
XA=BETA-DLOG v-DLOG ( 1 . D0-AA) ) /ALPHA 
RETURN 

4 BB=RAN ( I SEED ) 

SDX=DSQRT (DLOG < 1 . D0+XSTD**2) ) 

UX=DLOG (XMEAN) 

W=DSQRT (-2. D0*DLOG (AA) ) *T5COS (PI2*BB) *SDX+UX 
X=DEXP ( W) 

X A=DEXP (-Wh-2. D0*UX) 

RETURN 

5 X=BETA* (-DLOG ( AA) ) ** (-1 . D0/ ALPHA) 

XA=BETA* (-DLOG ( 1 . D0-AA) ) ** (-1 . DG/ ALPHA) 

RETURN 

END 

SUBROUTINE SECAl (COV, I SIGN, ALPHA) 

IMPLICIT DOUBLE PRECISION (A-H,0-Z> 

C I SIGN = 1 ; WE I BULL DIST. 

C = 2j FRECHET DIST. 

F (X , COV) =- ( 1 . D0+COV**2) *GAMMA (X) #*2+GAMMA (2. *X) 

IF ( ISIGN. EQ. 1 > X l=COV** ( 1 . 08) 

IF ( ISIGN. EQ. 2 ) XI =COV** ( . 677) /2. 33 
IF < ISIGN. EO. 2. AND. Xl.GT. .49D0) XI -.48999999 
7 IF ( ISIGN. EO. 1 ) F 1 -F (XI, COV ) 

IF ( ISIGN. EO, 2> F 1 =F < -X 1 , COV) 

IF (DABS <F1 ) .LE. 1 . D-10) GO TO 1 
X2=X 1+. 01DZ 

IF ( ISIGN. EO. 1 ) F2=F(X2,C0V) 

IF ( ISIGN. EO. 2) F2=F (-X2 , COV ) 

X X = X2 

10 CONTINUE 

IF (DABS (XX-X1) .GE. l.D-9) THEN 


Secant method for 
computing Weibull 
and Frechet exponents 


RAN is library 
uniform random 
number generator 
for CYBER 175 


GENX obtains random 
samples from the 
distribut ions 
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466: 
467: 
468: 
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470: 
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473: 
474: 
475: 


SUBROUT I NE E VD ( XME AN , STD , ALPHA , BETA , P I > 

IMPLICIT DOUBLE PRECISION <A-H,b-Z> 

ALPHA=PI / (STD*DSQRT (6. D0> > Computes EVD 

BETA=XMEAN-. 5772 1566490 153/ALPHA parameters 

RETURN 

END 


DOUBLE PRECISION FUNCTION GAMMA (Yl) 

IMPLICIT DOUBLE PRECISION <A-H,0-Z) 

COMMON /TWO/PI ,SPI2,PI2 1 

X*Y1 + 1 . D+0 8 annna 

Z = X function 

IF (X. GE. 6. 0D+0) GO TO 456 J __ 

N=INT (X) 

Z« (6. 0D+0) -N+X 
Y*1 . D+0/Z**2 

ALG* ( Z-. 5D+0) *DLOG ( Z ) +. 5D+0*DLOG (PI2>- 

Z- < 1 . D+0/ < 12. D+0*Z> )*( ( (Y/0. 14D+3-1.D+0/0. 105D+3)*Y+ 

1 . D+0/. 3D+2) *Y-1 . D+0) 

IF (X. GE. 6. D+0) GO TO 457 
ITE*6-N 
DO 3 J= 1 , ITE 
A=X*J-l.D+0 
ALG«ALG-DLOG (A) 

CONTINUE 
GAMMA=DEXP (ALG) 

RETURN 

END 


Note: The performance function must be introduced in subroutine LSFMC. 

For an example of subroutine LSFMC, see the last page of Appendix B. 
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Introduction 

The structural reliability analysis methods developed during the past 
ten to fifteen years can be used to establish the CDF of complicated 
structural response function by forming the so-called limit state function 
or performance function [1], In the application of these methods [2 - 4], 
the Taylor's series expansions are taken at the most probable points. For 
a given response function value, there Is a corresponding most probable 
point which needs to be found using proper optimization or Iteration 
algorithm. Because at each of the most probable point, there Is no error 
In the function and the error Is small around the most significant region 
for probability calculations, reasonably accurate solutions are assured. 
Indeed, experience has Indicated that the applications of these methods 
usually results in high quality CDF estimation. However, when the response 
function Is complicated, and the computations of the response variables are 
tedious, the above methods tend to be difficult to be Implemented and/or 
are prohibitively time-consuming. 

Presented here is a more efficient scheme Which Is suitable for 
constructing the cumulative distribution function (CDF) of any complicated 
function which has no explicit form, i.e., the objective function can not 
be expressed in algebraic form. The method is particularly suitable for 
the cases where the computation of the objective function Is time consuming 
such that the Monte Carlo method becomes prohibitively costly. 


Efficient Method of Constructing CDF using the Most Probable Points 

The efficient method of constructing the CDF of a function starts with 
a linear approximation of the response function about the mean values of 
the independent random variables. Then the CDF values and the associated 
most probable points for several "predicted" response function values will 
be computed. For any selected CDF value, however, the "predicted" response 
function is not accurate if the response function is nonlinear, therefore, 
the corresponding response function value will be "corrected" by solving 
the actual values at the predicted most probable point. 
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In order to show how the method works, a simple example Is established 
to detail the above procedure. The example Is a cantilever beam. The 
random variables Involved are the applied force, P, and the length, L, 
which are assumed to be normally distributed variables. The mean and 
standard deviation of P are (0.223, 0.019) lbs. and the mean and standard 
deviation of L are (20, 1) Inches. The maximum stress, S, at the fixed end 
of the beam Is: 


S = 787LP (1) 

The mean value of S Is approximately 3500 psl. 

Define the reduced variables u^ and u 2 as 

u x = (P - 0 . 223 ) /0 .019 (2) 

u 2 = (L - 20.)/l. (3) 


Thus 


P = 0.223 + 0.019u 1 ( 4 ) 

L = 20 + u 2 (5) 

By substituting equations 5 and 6 into equation 1, the stress becomes 


S = 3510 + 300uj + 175u 2 + 15ulu 2 


( 6 ) 


By assigning a value for S, an iso-stress curves can be plotted on a two 
dimensional u space as shown in Figure 1. Note that U} and u 2 are 
standardized normal variables (with zero mean and unity standard deviation) 
because P and L are normal variables. Therefore, on the uj,u 2 coordinate 
system, the joint probability density function Is rotational ly symmetric. 
The most probable point for a given S is easily identified as the point on 
a iso-stress curve which is nearest to the origin. 

Now we can start the approximation procedure by taking the first order 
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expansion of S about the mean. Usually, one win operate on the L, P space 
and then transform to uj, U 2 space. In this example, since the 
transformations are linear, we can use Eq. 6 directly. The first order 
expansion about the mean values (uj * 0 and u 2 = 0) results In 

S = 3510 + 300UJ + 175u 2 (7) 

This equation Is exact at u^ * u 2 = 0 (where S = 3510) only, but can be 
used as an approximation for other S values. 

Based on Eq. 7, S is also a normal variable with a mean of 3510 and a 
standard deviation of 347. It is obvious, however, that the accuracy of 
the CDF of S will depend on the truncated higher order terms. 

Traditionally, a low order expansion (such as eq. 7) Is only used to 
estimate the mean and the standard deviation. The CDF cannot be accurately 
approximated. 

For illustration purposes, only one CDF value will be considered. Let 
Sj = 3510 psl and S 2 = 4500 psl where Sj curve Ms linear and passes through 
the origin and S 2 is parallel to in the u space. S 2 may be called a 
"predicted" stress since eq. 7 is assumed to hold for all the stress values. 
The predicted S 2 curve has a most probable point which Is a point nearest 
to the origin. Assuming that eq. 7 is accurate, the first order 
reliability analysis method gives the following probability estimate: 

P(S > 4510 psi) = ®(- 8) (8) 

where 6 is the minimum distance. The approximation, however, Is not 
accurate because the most probable point derived was based on the 
Inaccurate S equation. In fact, by substituting the most probable point 
(derived by assuming S=4500 in Eq. 7) into Eq. 6, the exact value Is S = 
4660 psi. The iso-stress curves S^, S 2 and the exact S = 4660 are shown in 
Fig. 1. The exact curve is nonlinear and passing through the predicted 
most probable point. Since the predicted and the corrected curve match 
closely around the most probable region for S > 4660 (note that in Fig. 1, 
the two minimum distances are approximately equal.), the figure suggests 
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the following approximation: 

P(S Exact > 4660) * P($ Predicted > 4510) ( 9 ) 

Mathematical Formulations 

The above approximation can be formulated as follows. Let 

Z(X) * a 0 + ra i X i + E = Z p + E (10) 

where Z (X) is a function of the random variables, X. Z p is a random 
variable representing the sum of the Taylor's first order terms and E 
represents the sum of the higher order term. The error term should 
actually be a random variable, but in the present method It will be 
approximated by a deterministic value. E Is defined based on the most 
probable point, 1 .e. , 

E = Z(most probable point for Z p = z p ) - z p (11) 

where the most probable point Is defined as a set of values of X which 
maximize the joint probability density function of X subjected to the 
constraint that Z p (X) = z p . The most probable point can be found using the 
reliability analysis method [1]. 

Define the exact deterministic value of Z as z, then, 

P(Z > z) = P(Z p + E > z) 

= PU p > z - E) (12) 

- p ( z p > ^p) 


where z p is the predicted Z value using Z p . Equation 12 can be stated as : 
the probability of exceedance of Z p at z p is approximately equal to the 
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probability of exceedance of the exact Z at the value of 2 computed using 
the most probable point of Z p =z p . By replacing Z by S, z p by 4510, and z by 
4610, Eq. 12 becomes Eq. 9. 

To construct the entire CDF, the above demonstrated procedure can be 
repeated for other probability levels. Note that there is no limitation on 
the number of random variables and that the random variables can be any 
distribution. 

The first order Z p seems to be able to provide good approximation 
solutions as demonstrated In the following example. Improvements can be 
made by Including the second order terms In Z p . Alternatively, one can 
perform additional first order Z p analyses at the tall regions using the 
predicted most probable points. 


Establishing CDF - Example 

The above algorithm has been used successfully to establish a CDF of 
a problem. The problem is similar to the previous cantilever beam problem 
except that the thickness of the beam Is also modelled as a random 
variable. The goal is to estimate the CDF of the maximum stress. 

Figure 2 shows the resulting CDF based on the analytical solution of 
the stress. CDF curves are plotted on the normal probability paper (the 
CDF of a normally distributed variable will be a linear line on this 
paper); the Y coordinate uses u as the basic unit where u is a standardized 
normal variates. 

Using the conventional first-order-mean-expansion, the resulting CDF, 
in Fig. 2 , is nearly a straight line Indicating that S is approximately 
normal. This is because the approximating function Is linear and the 
random variables studied are normal or nearly normal. 

By applying the new algorithm, ten most probable points corresponding 
to ten CDF values are computed, and used to compute ten additional 
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deterministic solutions. These "new" stress values are the corrected 
values for the "old" CDF values. Figure 2 shows that the corrected CDF 
curve is very close to the "exact" (based on one million Monte Carlo 
simulations) CDF curve except at the region where u < -3, as shown In Fig. 
2 . 


In this example, the difference between Z e and Z p ranges from 0.2 % 
(for u = -0.5) of Z e to 32 % (for u * -4.3) of Z e , suggesting that the 
response function Is significantly nonlinear. This Is reflected by the 
fact that. In Fig. 2, the corrected CDF curve is significantly non-normal. 
Therefore, by using the new algorithm, It Is possible to assess the results 
by comparing Z e and Z p . Improvement is necessary only when the difference 
Is large. 

To improve the accuracy at the tail regions, there are two possible 
ways : (1). Take two more expansions at the tail regions (e.g., at u = 

-2.5 and u - 2.5). (2). Use quadratic or incomplete quadratic 

approximation about the mean values. The first method may be more 
appropriate when the quadratic approximations are difficult to obtain 
although the latter may provide more accurate results for problems 
involving highly nonlinear functions. 

The performance of the new algorithm has also been evaluated using 
Fig. 3. In this figure, the CDFs of the three mean-based approximations 
to the exact solution are constructed to compare with the exact CDF. The 
three approximations are: linear, "Incomplete" quadratic (second order 
mixed-terms are neglected), and (complete) quadratic expansions about the 
mean values of the independent random variables. By comparing the results 
of Fig. 3 with those of Fig. 2, it can be concluded that the new 
algorithm with only first order expansion performs better than the 
conventional quadratic expansion. Due to the fact that the complete 
quadratic approximations are much more difficult to obtain than the first 
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order approximations, the new algorithm with only first order approximation 
seems to be very suitable for estimating the CDFs for complicated 
functions. 


Summary 


The performance of the new algorithm using the demonstrated example Is 
excellent by noting that only a number of deterministic solutions. In 
additional to the first-order-mean-expansion, are required. The results 
suggest that the new procedure Is efficient and can be used to provide good 
CDF estimations for engineering analysis problems. 
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Figure 1. An Example showing the New CDF Estimation Algorithm 
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Figure 2. The Results of an Example using the New CDF Estimation Method 
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Figure 3. CDFs using Diffrent Approximation Levels 
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